Skip to main content


Differential expression of the aryl hydrocarbon receptor pathway associates with craniofacial polymorphism in sympatric Arctic charr

Article metrics



The developmental basis of craniofacial morphology hinges on interactions of numerous signalling systems. Extensive craniofacial variation in the polymorphic Arctic charr, a member of the salmonid family, from Lake Thingvallavatn (Iceland), offers opportunities to find and study such signalling pathways and their key regulators, thereby shedding light on the developmental pathways, and the genetics of trophic divergence.


To identify genes involved in the craniofacial differences between benthic and limnetic Arctic charr, we used transcriptome data from different morphs, spanning early development, together with data on craniofacial expression patterns and skeletogenesis in model vertebrate species. Out of 20 genes identified, 7 showed lower gene expression in benthic than in limnetic charr morphs. We had previously identified a conserved gene network involved in extracellular matrix (ECM) organization and skeletogenesis, showing higher expression in developing craniofacial elements of benthic than in limnetic Arctic charr morphs. The present study adds a second set of genes constituting an expanded gene network with strong, benthic–limnetic differential expression. To identify putative upstream regulators, we performed knowledge-based motif enrichment analyses on the regulatory sequences of the identified genes which yielded potential binding sites for a set of known transcription factors (TFs). Of the 8 TFs that we examined using qPCR, two (Ahr2b and Ap2) were found to be differentially expressed between benthic and limnetic charr. Expression analysis of several known AhR targets indicated higher activity of the AhR pathway during craniofacial development in benthic charr morphotypes.


These results suggest a key role of the aryl hydrocarbon receptor (AhR) pathway in the observed craniofacial differences between distinct charr morphotypes.


Adaptive diversification of craniofacial structures in teleost fish offers a remarkable opportunity to study niche specialization in vertebrates [1]. The formation of craniofacial elements begins at early stages of development through complex and dynamic molecular programs within the tissues bearing relevant patterning information [2]. At the forefront of efforts to understand the underlying developmental and evolutionary mechanisms, species such as zebrafish and cichlids have been extensively studied [3, 4]. Addressing similar questions in non-model species is becoming commonplace due to recent advances in molecular techniques such as large-scale gene expression (transcriptome) profiling and related bioinformatics [58]. Among salmonid fishes, the polymorphic Icelandic Arctic charr (Salvelinus alpinus) in Lake Thingvallavatn is an outstanding system to investigate intraspecific variation of craniofacial structures [9, 10]. Four distinct subpopulations of Arctic charr are readily distinguishable in this lake based on their genetic, ecological, behavioural and morphological differences [915]. Two benthic morphs, large benthivorous (LB) and small benthivorous (SB), have an overshot or subterminal mouth, whereas two limnetic morphs, planktivorous (PL) and piscivorous (PI), have a terminal mouth and a higher number of gill rakers (Fig. 1) [10]. Results obtained in laboratory rearing experiments have suggested a genetic basis as well as a combination of maternal and environmental factors contributing to this divergent trophic morphogenesis [1618]. We propose that these genetic differences should be detectable as differential transcriptional dynamics of genes involved in craniofacial skeletal development and be reflected later in different, morph-specific phenotypes [19].

Fig. 1

Drawings of adult heads of the four Arctic charr morphs in Thingvallavatn, Iceland. The two limnetic morphs, piscivorous (PI) and planktonivorous (PL) have terminal mouths and pointed snouts and the two benthic morphs, large benthivorous (LB) and small benthivorous (SB) have blunter snouts and more sub-terminal mouths. Scale bars represent 1 cm. (drawings by Eggert Pétursson)

In cichlids, transcriptional changes of major signalling pathways during early craniofacial development have been attributed to phenotypic novelties in trophic architecture [20, 21]. Associations between heterochronic shifts in the expression levels of individual skeletogenic genes and distinct craniofacial morphogenesis were indicated in several fish species [2224]. However, it is not always the case that expression differences of a specific gene result in changes in its protein levels and subsequent phenotypic effects [25]. Instead, looking for transcriptional differences in sets of functionally related genes with conserved co-expression connectivity might provide more reliable evidence concerning the involvement of genes in a specific biological process [26]. More recently, similar efforts have been made to define co-expression gene networks involved in adaptive phenotypic divergence of feeding structures in an East African cichlid species [27] and two salmonid species (Arctic charr and whitefish) [19, 28]. We previously identified a conserved co-expression network consisting of genes involved in extracellular matrix (ECM) organization and skeletogenesis, which was more expressed in the developing craniofacial elements of benthic than in limnetic Arctic charr morphs [19]. Many of these genes are known as downstream targets of several inter-connected regulatory pathways [19]. It is possible that the molecular mechanisms underlying induction of this gene network also affect (activate/repress) other genes and possibly co-expressed networks of genes involved in craniofacial development and skeletogenesis. The identification of groups of genes with tight positive or negative co-expression could provide clues about upstream molecular pathways [29, 30] and even causative genetic changes.

In this study, we investigated transcriptional profiles in developing heads of benthic and limnetic Arctic charr morphs focusing on the period when key elements of the feeding structures are being formed, and, notably, a period when the basis of the benthic–limnetic distinction may be laid [19]. Using transcriptome data of contrasting Arctic charr morphs and through a stepwise approach assisted by co-expression data from mammalian species, we identified genes with higher expression during the morphogenesis of limnetic charr heads. Furthermore, to identify potential upstream regulators, we performed knowledge-based motif enrichment analyses on the regulatory sequences of a network of co-expressed genes (including the identified genes) to predict binding sites for known transcription factors (TFs). Finally, we found Ahr2b and other components of a signalling pathway mediated by aryl hydrocarbon receptor (AhR) to be differentially expressed between benthic and limnetic charr, which is consistent with the AhR pathway being a regulator of morph-specific, transcriptional dynamics. These results provide a glimpse into the regulatory networks that seem to shape craniofacial traits in Arctic charr and possibly other teleosts.


Fish stocks and sampling

We sampled sexually ripe adults of the four sympatric morphs from Lake Thingvallavatn; planktivorous (small limnetic) (PL), piscivorous (PI), small benthivorous (SB) and large benthivorous (LB) morphs, and a limnetic-like aquaculture stock (AC) from the Hólar breeding programme [31]. For each morph, eggs from several females were pooled and fertilized using milt from several males. Eggs were reared at approximately 5 °C in hatching trays (EWOS, Norway) under constant water flow and in complete darkness at Holar College experimental facilities in Verið, Sauðárkrókur. The water temperature was recorded twice daily and the average was used to estimate the relative age of the embryos using tau-somite (τ s) units, defined as the time it takes for one somite pair to form at a given temperature [32]. All sampling from the wild as well as rearing was performed according to Icelandic law and with appropriate permission.

Genes, databases and overrepresentation analysis

Based on an extensive literature survey and analyses of transcriptome data from early embryonic stages of contrasting Arctic charr morphs (unpublished), we selected our initial sets of candidate genes (see “Results” section) that play a role in craniofacial morphogenesis and/or skeletogenesis, for gene expression profiling. We also checked whether the orthologues of the candidate genes show co-expression relationship in data from mammalian species using COXPRESdb ( version 6.0 [33]). To predict the potential regulators of the genes with reduced expression levels during benthic head development, TF enrichment analysis was conducted using the list of >200 genes co-expressed with Chd4, Cspp1 and Dlg1 in humans and mice, as well as the WEB-based GEne SeT AnaLysis Toolkit (WebGestalt) v2 [34]. The involvement of the AhR was investigated by profiling the expression of 12 genes responsive to the pathway which also have been shown to be expressed during craniofacial development in fish and/or involved in vertebrate skeletal formation (Additional file 1: Table S1).

RNA isolation, cDNA synthesis and primer design for qPCR

We studied embryos of five Arctic charr morphs, collected at 6 tau-somite stages (178, 200, 216, 238, 256 and 275 τ s). This is the period during which the morphs start developing their distinct craniofacial morphologies [35]. The time points span early craniofacial cartilage formation, 178 (τ s), to a pre-hatching stage, when ossification of key trophic elements of the jaws has started at 275 (τ s). The embryos were preserved in RNA-later solution (Ambion) at −20 °C and later dechorionated and decapitated by applying a scalpel vertically in front of the pectoral fin under the light microscope (Leica S6E). For each morph and time-point, two RNA extraction replicates were obtained, each one deriving from six heads. Extractions were conducted in TRI Reagent (Sigma). The heads were homogenized with a disposable Kontes Pellet Pestle Cordless Motor tissue grinder (Kimble Kontes). RNA was extracted according to the manufacturer’s instructions, and DNA contamination was removed by DNase treatment (New England Biolabs). The quantification and quality control of RNA were performed using a NanoDrop ND-1000 UV/Vis-Spectrophotometer (NanoDrop Technologies) and agarose gel electrophoresis. cDNA was synthesized with 1 µg of RNA in a total reaction volume of 20 µl using the High Capacity cDNA RT kit (Applied Biosystems). cDNA was diluted threefold in nuclease-free water for further use in qPCR and several samples without addition of reverse transcriptase were prepared to confirm the absence of genomic DNA. An assembly of the Arctic charr transcriptome [36] was used for qPCR primer design. It is important to note that a whole genome duplication event occurred before the radiation of the salmonid family producing paralogues, most of which are retained in Arctic charr [37]. In this study, we were not focusing on differences between paralogues in expression or function and therefore designed primers to be non-paralogue specific. The only exception is Ahr2b where we looked at differential expression of two paralogues as discussed below. To identify the conserved exon/intron boundaries of the candidate genes, the assembled charr contigs were aligned to the genomic sequences of salmon orthologs (retrieved from the salmon database, SalmonDB [38]) using the NCBI Spidey software ( Primers were designed using Primer Express 3.0 software (Applied Biosystems, Foster City, CA, USA) and checked for self-annealing, hetero-dimers and hairpin structures by OligoAnalyzer 3.1 (Integrated DNA Technology) (Additional file 2: Table S2).

Real-time quantitative PCR and analysis of expression data

Real-time PCR was performed in 96-well PCR plates on an ABI 7500 real-time PCR System (Applied Biosystems) using Fermentas Maxima SYBR Green qPCR Master Mix (2x) (Fermentas). Each biological replicate was run in duplicate together with no-template control (NTC) in each run for each gene. The PCR was run with a 2 min hold at 50 °C and a 10 min hot start at 95 °C followed by the amplification step for 40 cycles of 15 s denaturation at 95 °C and 1 min annealing/extension at 60 °C. A dissociation step (60–95 °C) was performed at the end of the amplification phase to identify a single, specific melting temperature for each primer set. Primer efficiency values (E) were calculated through the LinRegPCR v11.0 programme ( [39] analysing the fluorescence data from the exponential phase of PCR amplification for each primer pair (Additional file 2: Table S2) and the cut-off was 0.9. The difference between Cq values (ΔCq) of the reference genes and the target genes was calculated for each gene t; ΔCqtarget = Cqtarget − Cqreference. The geometric mean of Cq values of two validated craniofacial reference genes, Actb and If5a1, was used for ΔCq calculations [40]. All samples were then normalized to the ΔCq value of a calibrator sample to obtain a ΔΔCq value (ΔCqtarget − ΔCqcalibrator). For visual comparisons of expression levels among the developmental time-points and morphs, the biological replicate with the lowest expression (highest ΔCq) was used as a calibrator sample. Relative expression quantities (RQ) were calculated based on the expression level of the calibrator sample (E −ΔΔCq) [41]. The RQ values were then transformed to logarithmic base 2 values (or fold differences; FD) [42] for statistical analysis. A two-way ANOVA followed by Tukey’s HSD post hoc tests were implemented to examine the effects of morph, time, and morph-by-time interaction on expression of the candidate genes as well as morph-specific expression difference for each gene. To assess similarity in expression patterns of the genes, Pearson correlation coefficients (r) were calculated for all gene pairs using the data from 5 morphs at 6 time-points (df = 28). To study the relationships between putative regulators and candidate genes, we calculated the correlation (Pearson and Kendall) between all candidate genes (including AhR targets) and the 9 potential regulators (8 predicted TFs and Foxq1, an AhR-regulated TF). Next, we tested whether the correlation coefficients varied between TFs or gene categories (higher, lower or equally expressed in benthic morphs), using ANOVAs and non-parametric. R ( was used for all statistical analysis [43].

Whole-mount in situ hybridization (WISH)

Whole-mount in situ hybridization was performed following a standard procedure for Atlantic salmon [44]. Embryos from a time-point (238 τ s) representing the early craniofacial bone and cartilage formation were fixed in 4 % (m/v) paraformaldehyde/PBS, dehydrated in a graded methanol series and stored in 100 % methanol. Primers designed for cDNA of selected Arctic charr genes (Ahr2b, Apaf1, Cfl1, Foxq1, Jup, Chd4, Cspp1, and Dlg1) generated PCR products of around 0.4–0.6 kbp (Additional file 2: Table S2) which were cloned into pCR4-TOPO vector (Invitrogen) and transcribed to antisense/sense digoxigenin (DIG)-labelled cRNA probes with T3 or T7 polymerases (Roche). Five or six embryos from each of the two morphs (PL and SB) were used for in situ hybridization at every relative age. After rehydration and dechlorination, embryos were treated with 20–40 µg/ml proteinase K (New England Biolabs) for 22–60 min, depending on the relative ages. The hybridization was performed with 1 µg/ml (DIG)-labelled RNA probes at 70 °C for 12 h, after which the embryos were incubated with alkaline phosphatase-conjugated anti-digoxigenin antibody (Roche) at 4 °C overnight. The hybridization signals were visualized after 3–12 h incubation at 4 °C using NBT/BCIP (Roche), depending on the expression levels of the genes. The specificity of the antisense probes was also verified by running control experiments with sense probes. Samples were imaged on a Leica MZ10 binocular microscope (Leica).


A group of transcriptionally correlated genes shows higher expression in the head of limnetic than benthic morphs

We began by studying the expression of six genes (Brg1, Carm1, Chd4, Kat6a, Med12 and Rbp2), which are known to play a role in early craniofacial morphogenesis and/or skeletal differentiation [4551]. Preliminary mRNA-seq data indicated higher expression of these genes during early embryonic stages in limnetic than in benthic morphs (Gudbrandsson et al. unpublished), and according to COXPRESdb [33] they show co-expression in mammalian species (Fig. 2a). By further examining the co-expression of these genes with other genes from mammals, we found six more genes (Cbp, Chd7, Cpsf1, Eftud2, Hcfc1 and Whsc1), which show conserved co-expression and are all known to play a role in craniofacial skeletal morphogenesis [46, 5256] (Fig. 2a). To characterize the expression of these 12 candidate genes in the developing head of Arctic charr, we profiled their relative expression in the heads of five morphs (SB, LB, PI, PL and AC) with qPCR, at six consecutive time points, spanning early craniofacial bone and cartilage formation to the prehatching stage (Fig. 2b). The expression of all genes, except Cbp, varied significantly over time, showing a general trend of gradual reduction during head development. All genes except Cpsf1 and Eftud2 were also expressed at different levels between morphs. However, only three genes, i.e. Chd4, Hcfc1 and Whsc1, showed lower expression in both benthic morphs (SB and LB) compared to the limnetic morphs (AC, PI, PL) (Fig. 2c). The benthic–limnetic expression difference was markedly larger for Chd4 than the two other genes within several developmental time-points (Fig. 2b). A positive correlation of expression levels between most of the candidates, including Chd4, Hcfc1 and Whsc1, was verified in Arctic charr (Fig. 2d).

Fig. 2

Expression of Chd4, Hcfc1 and Whsc1 was higher in developing heads of limnetic Arctic charr morphs. a Expression correlation of the first twelve candidate genes based on mutual rank (MR) analysis from COXPRESdb in mammalian species. The stronger correlations represented in darker green shading and with one, two and three asterisks indicating <2000, <1000 and <500 MR numbers, respectively. b Relative expression of the twelve genes in developing heads of SB, LB, PI, PL and AC at six stages. Error bars represent standard deviation calculated from two biological replicates where each biological replicate contains a homogenate of six heads. c Relative expression ratios were subjected to an analysis of variance (ANOVA) to test for expression differences among five Arctic charr morphs and six stages (M morph; T stage (time); MxT morphs by time effects; P values of <0.05, 0.01 and 0.001 are indicated by one, two and three asterisks, respectively). Subsequently, a post hoc Tukey’s HSD test was performed to analyse the expression of candidates among the morphs. Green to red colour gradient of morphs represent low to high expression levels and morphs with no connecting line have significantly different expression (α = 0.05). The bold and underlined genes displayed lower expression levels in the benthic morphs. d Pearson correlation coefficient (r) was used to assess the pairwise expression similarity between the candidate genes during craniofacial development. Blue represents positive expression correlation (the colour gradient showing critical values of r, 2-tail; df = 28). P values of <0.05, 0.01 and 0.001 are indicated by one, two and three asterisks, respectively

With the aim of identifying more co-expressed genes with potentially reduced expression in the heads of benthic morphs, we extended the expression analysis by including eight more genes. Seven genes were selected based on co-expression with Chd4, Hcfc1 and Whsc1 in mammals (using COXPRESdb [33]) (Fig. 3a) and involvement in vertebrate craniofacial development [5763]. An eighth candidate, Rev3l, was added to the study due to its vital role in early embryonic patterning and craniofacial formation [64] and co-expression with four of the new candidates in mammals (although not Chd4, Hcfc1 and Whsc1) (Fig. 3a). We found that expression of all eight genes differed significantly between morphs and over time, and four of them, i.e. Cspp1, Dlg1, Kdm5c and Rev3l, showed lower expression in both benthic morphs compared to all the limnetic morphs (Fig. 3b, c). The benthic–limnetic expression difference was larger for Cspp1 and Dlg1 than the two other genes within several developmental time points (Fig. 3b). Again the expression of most of these eight genes was positively correlated with Chd4, Hcfc1 and Whsc1 from the previous step (Fig. 3d). Taken together, these results show significantly higher expression of 7 genes out of the 20 tested in developing heads of limnetic morphotypes as well as their positive expression correlation during the developmental period under study.

Fig. 3

Expression of Cspp1, Dlg1, Kdm5c and Rev3 l was higher in developing heads of limnetic than benthic Arctic charr morphs. a Expression correlation of the additional eight candidate genes based on mutual rank (MR) analysis from COXPRESdb in mammalian species. The stronger correlations represented in darker green shading and with one, two and three asterisks indicating <2000, <1000 and <500 MR numbers, respectively. b Relative expression of the genes in developing heads of SB, LB, PI, PL and AC at six stages. Error bars represent standard deviation calculated from two biological replicates where each biological replicate contains a homogenate of six heads. c Relative expression ratios were subjected to an analysis of variance (ANOVA) to test for expression differences among five Arctic charr morphs and six stages (M morph; T stage (time); MxT morphs by time effects; P values of <0.05, 0.01 and 0.001 are indicated by one, two and three asterisks, respectively). Subsequently, a post hoc Tukey’s HSD test was performed to analyse the expression of candidates among the morphs. Green to red colour gradient of morphs represent low to high expression levels and morphs with no connecting line have significantly different expression (α = 0.05). The bold and underlined genes displayed lower expression levels in the benthic morphs. d Pearson correlation coefficient (r) was used to assess the pairwise expression similarity between the candidate genes during craniofacial development. Blue represents positive expression correlation (the colour gradient showing critical values of r, 2-tail; df = 28). P values of <0.05, 0.01 and 0.001 are indicated by one, two and three asterisks, respectively

Potential upstream regulators of the co-expressed genes

With the objective of finding potential upstream regulatory pathway(s) and/or TFs influencing the co-expressed genes, we selected the three genes, Chd4, Cspp1 and Dlg1, with persistent and strong reduction in expression in developing heads of benthic morphs and with positive co-expression in mammals and Arctic charr. We retrieved a list of over 200 genes showing co-expression in mammals with these three genes (Additional file 3: Table S3). This list was used as input for a knowledge-based TF enrichment analysis in human and mouse using WebGestalt v2 [34]. In both human and mouse, binding sites belonging to 17 TFs were significantly overrepresented on promoters of the input genes (Additional file 4: Table S4). We chose the eight most significantly enriched TFs for further gene expression analysis (Fig. 4a). The results revealed differential expression of all the eight TFs over time and seven TFs between the morphs (Fig. 4b). However, only two TFs (Ahr2b and Ap2) were differentially expressed between benthic and limnetic charr, with higher levels of expression for Ahr2b in benthic morphs and higher expression of Ap2 in limnetic morphs. This implicates Ahr2b and Ap2 as potential transcriptional regulators (i.e. transcriptional repressor and activator, respectively) of the co-expressed genes identified above.

Fig. 4

Differential expression of Ahr2b and Ap2 in developing heads of Arctic charr morphs. a Relative expression of eight selected TF genes in developing heads of SB, LB, PI, PL and AC at six stages. Error bars represent standard deviation calculated from two biological replicates where each biological replicate contains a homogenate of six heads. b Relative expression ratios were subjected to an analysis of variance (ANOVA) to test for expression differences among five Arctic charr morphs and six stages (M morph; T stage (time); MxT morphs by time effects; P values of <0.05, 0.01 and 0.001 are indicated by one, two and three asterisks, respectively). Subsequently, a post hoc Tukey’s HSD test was performed to analyse the expression of candidates among the morphs. Green to red colour gradient of morphs represents low to high expression levels and morphs with no connecting line have significantly different expression (α = 0.05). The bold and underlined gene(s) displayed distinct benthic–limnetic expression dynamics

Differentially expressed AhR target genes show transcriptional correlation with the previously identified genes

The observations above led us to hypothesize that the signalling pathway mediated by aryl hydrocarbon receptors (AhR pathway) is involved in the emergence of distinct benthic–limnetic craniofacial morphologies. This hypothesis is also bolstered by the fact that a homozygous loss-of-function but viable mutant of an aryl hydrocarbon receptor (Ahr2) in zebrafish has a pointed snout and protruding lower jaw [65]. The phenotypic differences between the mutant and a wild-type zebrafish to some extent reflect the phenotypic dichotomy seen in Arctic charr, i.e. the long lower jaw and pointed snout of limnetic morphs versus the shorter lower jaw and blunt snout of benthic morphs [65]. In view of this, we selected twelve genes for expression analysis that are known downstream targets of the AhR pathway and also possibly involved in vertebrate craniofacial morphogenesis and/or skeletogenesis (Additional file 1: Table S1). Interestingly, the qPCR results indicated that eight of these genes had increased expression (and one gene had decreased expression) in the heads of benthic morphs (relative to limnetic morphs) during development (Fig. 5). The magnitude of the expression differences varied considerably among the genes with Cyp1a1 and Foxq1 showing the largest expression differences between morphs in several developmental time-points (Fig. 5a). We calculated Pearson’s correlation coefficients of the expression levels for AhR target genes showing distinct benthic–limnetic transcriptional dynamics together with the previously identified genes, over all morphs and time-points (Fig. 6a) and found strong positive expression correlation between pairs of AhR targets, and also between them and Ahr2b. On the other hand, the same genes showed negative expression correlation with genes identified in previous steps (mostly with Chd4, Cspp1 and Dlg1) (Fig. 6a). Strikingly, the only AhR target gene with reduced expression in benthic morphs, Apaf1, displayed positive correlation in expression with all the genes found in former steps as well as with the Ap2 gene (Fig. 6a). Ap2 showed positive expression correlations with the set of genes that had higher expression in limnetic heads, but negative expression correlations with Scin and Foxq1, which is itself an AhR-regulated transcription factor.

Fig. 5

Differential expression of AhR-regulated genes in developing heads of Arctic charr morphs. a Relative expression of twelve AhR-regulated candidate genes in developing heads of SB, LB, PI, PL and AC at six stages. Error bars represent standard deviation calculated from two biological replicates where each biological replicate contains a homogenate of six heads. b Relative expression ratios were subjected to an analysis of variance (ANOVA) to test the expression differences amongst five Arctic charr morphs and six stages (M morph; T stage (time); MxT morphs by time effects; P values of <0.05, 0.01 and 0.001 are indicated by one, two and three asterisks, respectively). Subsequently, a post hoc Tukey’s HSD test was performed to analyse the expression of candidates among the morphs. Green to red colour gradient of morphs represents low to high expression levels and morphs with no connecting line have significantly different expression (α = 0.05). The bold and underlined gene(s) displayed differential benthic–limnetic expression dynamics

Fig. 6

Significant expression correlations between Ahr2b, AhR-regulated genes and the seven previously identified genes and their distinct benthic–limnetic expression dynamics. a Correlation analyses revealed negative or positive expression correlation of Ahr2, AhR-regulated genes and the seven previously identified genes during craniofacial development. Pearson correlation coefficient (r) was used to assess the pairwise expression similarity between the candidate genes. Blue shading represents positive and red shading represents negative expression correlation and the colour gradients indicate correlation coefficients (r) above critical values (2-tail; df = 28). P values of <0.05, 0.01 and 0.001 are indicated by one, two and three asterisks, respectively. b Analyses of correlations (Kendall’s tau) between the expression of Ahr2, Ap2 and Foxq1 and the expression of other genes classified in three groups: an Up group including genes with significantly higher expression in the benthic charr morphs, a Down group with genes with significantly lower expression in benthic morphs and a Rest group including genes without significant expression difference between benthic and limnetic charr. Non-parametric Kruskal–Wallis ANOVA indicated difference among gene groups for the three TFs (P < 0.005). c A distinct morph and time-dependent expression pattern of the AhR-regulated and the seven previously identified genes were seen in developing heads of Arctic charr. The selected genes were subjected to double clustering based on expression correlation of the genes and between morphs by time effects across six stages. Both gene and morph-time clustering dendrograms include two main branches and each branch contains a few sub-clusters. A clear benthic–limnetic separation is observed in sub-clusters of the morph-time dendrogram (benthic samples depicted in red). Also, a clear separation of AhR-induced and the seven identified genes is observed in the gene clustering. Blue and yellow shadings represent higher and lower relative expression, respectively

To further assess expression associations between the TFs and the other genes studied here, we calculated the correlation between expression levels of all the genes together and the TFs (Fig. 6b, Additional file 5: Figure S1). The genes were classified into three categories, “Up”, “Down” and “Rest” based on significant expression differences between benthic and limnetic morphotypes. Genes in the Up group had higher expression in benthic charr, while the Down group showed lower expression in benthic charr. The Rest category consists of genes with no significant expression differences between the two morphotypes. As was predicted, Ahr2b showed significant positive correlation with its target genes (Up-genes in Fig. 6b) and a negative, but weak relationship with the Down-genes group. Ahr2b levels did not associate on average with genes in the Rest group. The Down genes were positively correlated with six of the other candidate regulators (Ap2, E2f1, Maz, Nrf1, Sp1 and Yy1) (Fig. 6b, Additional file 5: Figure S1), of which only Ap2 differed significantly between the morphotypes. The opposite effect, i.e. negative correlation with the Up genes, was generally not observed for these other TFs. Curiously, however, these factors all showed positive correlation with genes of the Rest category. The correlation of Foxq1 expression with the three groups of genes was similar to that seen for Ahr2b, except Foxq1, which is a TF known to be activated by AhR signalling, showed stronger negative association with the Down genes (Fig. 6b).

To summarize the overall patterns of expression, we conducted double hierarchical clustering of relative expression values of genes vs. morph and time (Fig. 6c). The genes are divided into two main branches. The AhR targets and Ahr2b (in bold) clustered in one branch and the rest of the genes with reduced expression in benthic morphs (including Ap2 and Apaf1) clustered in the other branch. Predictably, the clustering of samples (morph-time) revealed a separation of morphotypes with benthic morphs (red bar) on one of the two main branches and the limnetic morphs (black bar) on the other. The exceptions to this were PL and PI at early time-points, and AC which clustered with the benthic samples at the last time-point (Fig. 6c). Taken together, these results suggest higher activity in the AhR signalling pathway, probably through Ahr2b, during the craniofacial development of benthic Arctic charr morphs.

Craniofacial expression patterns of differentially expressed genes

To determine the spatial craniofacial expression and differences in expression pattern between benthic and limnetic morphs, we investigated in situ expression of eight genes. In addition to the differentially expressed Ahr2b, we selected the three AhR-regulated genes with higher expression level in the benthic morphs (Cfl1, Foxq1 and Jup) and four genes with reduced expression in the benthic morphs (Apaf1, Chd4, Cspp1 and Dlg1). The spatial expression pattern of these genes was studied in two contrasting morphs (PL and SB) at time point 238 τ s during early craniofacial skeletal morphogenesis. In concordance with the correlation analyses, all seven genes showed craniofacial expression (Fig. 7). A similar expression pattern could be observed for the selected genes in developing heads of both PL and SB. Five of the genes, Ahr2b, Apaf1, Cfl1, Foxq1 and Jup, had pronounced expression in regions surrounding the mouth, ventral and ventrolateral facial elements and pharyngeal arches. In the case of Cfl1 and Foxq1, which showed large expression differences between the morphs in the qPCR analyses (Fig. 5), a more pronounced and extended expression was detectable in ventral facial elements and pharyngeal arches of the SB morph (marked area in Fig. 7). The most evident spatial differences in expression were observed for Cfl1, where in the SB morph the expression was extended all over the region encompassing frontonasal region, the lower jaw, ventral facial elements and pharyngeal arches (Fig. 7). The expression of Chd4, Cspp1 and Dlg1 did not show clear spatial patterns, but was diffused all over the head mesenchyme without any noticeable spatial differences between the morphs.

Fig. 7

Craniofacial expression patterns of candidate genes in developing heads of two contrasting Arctic charr morphs (SB and PL). In situ hybridization shows the overlapping expression pattern of Ahr2b, Apaf1, Cfl1, Foxq1 and Jup at the 238 (τ s) time-point. Ventral and lateral views indicate the overlapping expression of the genes in the area surrounding the mouth, lower jaw and pharyngeal arches. An extended expression pattern is observed for Cfl1 and Foxq1 in ventral facial elements and pharyngeal arches of SB morph (marked area). Furthermore, Cfl1 had higher expression in the fronto-nasal region in SB (marked area). Chd4, Cspp1 and Dlg1 showed diffuse expression in most parts of the head mesenchyme without noticeable spatial expression differences between the two morphs


An increasing number of gene expression studies have been launched to unravel the molecular mechanisms underlying adaptive morphological novelties in feeding structures of a variety of fish species [1923, 27, 28, 66, 67]. Developmental shifts in the relative timing of skeletal gene expression, i.e. collagen genes, have been shown to underpin the radiation of feeding structures in Antarctic notothenioid fish along the benthic–limnetic axis of habitat [22]. Similarly, a heterochronic shift in developmental expression of a Calmodulin gene has been suggested to drive the variable growth rates of the jaws across Belonoid fishes [23]. In addition to the heterochrony in expression, the developmental differences in the level of gene expression have been repeatedly linked to occurrence of adaptive skeletal changes in teleost fish [1921, 28, 68, 69]. It is important to note that a highly similar skeletal phenotype might emerge from the altered expression of different regulatory genes as two distinct transcriptional mechanisms give rise to the loss of pelvic skeleton in stickleback and fugu [69, 70]. Studies in cichlids have shown that genetic changes in regulatory elements may account for altered expression of genes involved in skeletogenesis and the subsequent morphological variations [20, 71, 72]. Also, external stimulants during early development, such as mechanical strain in bone, can alter the expression of skeletal genes (e.g. ECM genes) and result in adaptive phenotypic plasticity in feeding structures [27, 66]. The molecular basis of adaptive morphological novelties in feeding structures can be further elucidated through the expression analysis of functionally and transcriptionally related gene networks. The identification of such gene network provides a promising foundation to investigate potential upstream regulatory elements and interacting signalling pathways [26, 29, 30]. The related examples include the identification of gene expression networks associated with evolution of oral jaw dentition [73] and pharyngeal jaw developmental plasticity [27] in East African cichlids, as well as craniofacial skeletal diversity in sympatric Arctic charr [19] and lake whitefish [28].

Recently, we identified a gene network involved in ECM organization, cell-substrate adhesion and skeletogenesis with higher levels of expression in developing heads of benthivorous Arctic charr than in limnetic morphs [19]. In the same study, we predicted several transcription factors, such as Ap1 and Ets2, as potential regulators of the network. Although these TFs could be key elements causing the induction of the ECM-related gene network, they are likely to be regulated through pathways further upstream. Such upstream pathways might also regulate other gene networks involved in craniofacial morphogenesis. In general, the activation of a particular signalling pathway leads to both induction and repression of different sets of downstream genes during embryonic development [29]. Therefore, we set out to search for co-expressed genes involved in skeletogenesis/craniofacial formation focusing on genes that had higher expression during head development in limnetic than in benthic morphs of Arctic charr. Such genes may be a part of networks that are transcriptionally linked to the previously described ECM-related gene network through common upstream regulatory elements.

In an iterative selection for co-expressed genes, we found seven candidate genes, Chd4, Cspp1, Dlg1, Hcfc1, Kdm5c, Rev3 l and Whsc1, showing reduced expression in benthic morphs. In a recent study on mammalian embryonic cells, Chd4 has been shown to be a potent co-repressor of sry-box containing gene 9 (Sox9), a highly conserved gene involved in skeletal development [47]. While the role of Chd4 in craniofacial development has not been investigated, other genes such as Cspp1, Dlg1, Hcfc1 and Whsc1 are well known to be essential for normal development of craniofacial skeletal elements [5557, 59]. In our study, the positive expression correlation of these genes could suggest their co-regulation through shared transcriptional regulatory elements during Arctic charr head development. Using co-expression data from mammals [33], we found over 200 genes to be co-expressed with Chd4, Cspp1 and Dlg1 as well as with at least one of the other identified genes. Binding sites for 17 TFs were significantly enriched on the promoter regions of these genes in human and mouse. By profiling the expression of eight predicted TFs with most significantly overrepresented binding sites, we found two of them, Ahr2b and Ap2, to be differentially expressed between benthic and limnetic morphotypes. Curiously, a zebrafish Ahr2 mutant displays an extension of the ethmoid and mandibular regions, shallower head and longer snout which is a reminiscent phenotype of limnetic Arctic charr morphology [16, 65]. Furthermore, Ap2-dependent signals from epithelium have been demonstrated to induce skeletogenesis during craniofacial development in zebrafish [74]. It is unclear, however, whether Ahr2 and Ap2 show a direct transcriptional interaction during development.

Aryl hydrocarbon receptor (Ahr) belongs to the basic-helix-loop-helix (bHLH)/Per-Arnt-Sim (PAS) family of heterodimeric transcriptional regulators and it was originally described as a receptor that mediates many of the toxic effects of environmental contaminants of the dioxin family [75]. It is highly unlikely that the actual role of AhR is to sense environmental contaminants, but the search for endogenous ligands is still ongoing (reviewed in [76]). There is, however, increasing evidence that the AhR pathway plays important roles in a variety of developmental, physiological and immunological processes [75, 7781]. The pathological effects of AhR pathway induction in zebrafish are exerted mostly through the transcriptional regulation of a large number of downstream genes [82]. It is also known that the activation of the AhR pathway not only alters the expression of ECM remodelling genes, but the pathway itself can also be affected by changes in ECM organization and cell-substrate adhesion [83]. Many of the identified TF binding sites in this study (such as binding sites for Ahr, Ap2, E2f1, Max, Myc, Myb, Nf1, Sp1, Usf, and Yy1) (Additional file 4: Table S4) are also predicted on the promoter of AhR-regulated genes in mammals [78, 84]. This suggests that these TFs cooperate in the regulation of AhR target genes that could be conserved across vertebrates. To further investigate the involvement of the AhR pathway, we examined the transcriptional dynamics of 12 downstream targets of the pathway which were selected based on existing literature of their craniofacial expression pattern and/or their role in vertebrate skeletal formation (Additional file 1: Table S1). All of these candidates, except Apaf1, are known to be induced by activation of the AhR pathway. Accordingly, we found higher expression of most candidates in developing heads of benthic morphs, whereas only Apaf1 displayed reduced expression relative to limnetic morphs (known to be down regulated by AhR). Most of the AhR-regulated candidates as well as Ahr2b had negative expression correlation with the seven genes identified at the beginning of this study and only Apaf1 showed positive expression correlations with them (Fig. 6a). The AhR-regulated candidates also had positive expression correlation with the previously identified co-expression network of ECM genes during the same developmental time-points ([19], data not shown). These results suggest higher activity of the AhR pathway during head development of benthic than limnetic morphs.

In Atlantic salmon, multiple aryl hydrocarbon receptor genes are reported, four of which are paralogues of Ahr2, with markedly higher expression than the two paralogues of Ahr1 [85]. Here, we only report results for one Ahr2 paralogue (Ahr2b) which had the highest expression level in a developmental transcriptome of Arctic charr [36]. Two other paralogues (Ahr2a and Ahr2c) were not detected in the transcriptome data, whereas the third one (Ahr2d) was identified with lower expression but qPCR analysis with paralogue specific primers showed no differences in expression between benthic and limnetic morphs (data not shown). Developmental effects of AhR pathway activation on skeletogenesis are mediated by Ahr2 in zebrafish and transcriptional suppression of Ahr2 could partially rescue the resulting inhibition of jaw chondrogenesis during development [86]. Analysis of the non-functional Ahr2 mutant in zebrafish illustrated the endogenous role of the AhR pathway in craniofacial development [65], but the molecular mechanism for this craniofacial effect has remained unclear. So far, investigation of potential downstream targets in zebrafish identified a contribution of Sox9 to AhR-induced jaw malformation [87]. In a previous study, we did not detect such differential expression of Sox9 at early stages of craniofacial skeletal formation in different Arctic charr morphs [40]. In zebrafish, down-regulation of the Sox9 isoform (Sox9b) does not occur immediately after induction of AhR pathway which might suggest an indirect mechanism for AhR regulation of Sox9b [87]. It has been reported that Ap2 negatively regulates chondrocyte differentiation by suppressing Sox5/6, but not Sox9 [88]. In contrast to Ahr2b, the higher expression of Ap2 in limnetic morphs and the reduced jaw skeleton and shorter snout of Ap2 mutant in zebrafish [65, 74] suggest counteracting roles for Ap2 and Ahr2b in mediating differential effects on craniofacial skeletogenesis.

Cyp1a1 is known as the most potently induced AhR target [75], and its expression has been identified in craniofacial tissues of teleost fish species [89]. In our study, Cyp1a1 shows strong benthic–limnetic expression differences (Fig. 5). A functional study revealed that suppression of Cyp1a1 expression does not prevent craniofacial malformation induced by activation of AhR during zebrafish development [90]. This demonstrates that even though induction of Cyp1a1 is an indicator of AhR pathway activation, it is not the key target of the pathway mediating the effects on craniofacial skeletal formation.

Among the other AhR downstream targets which could be involved in mediating the craniofacial effects (Additional file 1: Table S1), we found 8 more genes showing benthic–limnetic differential expression over the examined developmental period (Fig. 5). Three of the genes, Apaf1, Cfl1 and Foxq1, might be of importance for the development of the craniofacial characteristics induced by AhR pathway activity. Apaf1 is an evolutionarily conserved cytosolic protein with a critical role in regulation of developmental apoptosis in mammals [91] and we found that its expression is higher in developing heads of limnetic morphs. A non-functional Apaf1 mutant in mice displays defective craniofacial development through altered regulation of hedgehog (Hh) signalling [92]. The Hh pathway itself is known to be a key driver for adaptive variation of craniofacial skeletal structures in cichlids [93]. In mammalian cells, the expression of Apaf1 was shown to be inhibited by overexpression of Ahr through a complex formed by AhR and E2f1 which binds to the Apaf1 promoter at a region containing E2f1 binding site, but no AhR binding sites [94]. In general, AhR can act as a co-repressor of E2f-dependent transcription, and this is particularly interesting because binding sites for E2f were significantly enriched in our TF overrepresentation analysis [95].

Cfl1 encodes a protein (cofilin-1, non-muscle) that plays a role in the regulation of cell morphology and cytoskeletal organization by affecting the polymerization of actin [96]. Activation of the AhR pathway can induce the expression of Cfl1 [97]. Moreover, loss of function mutations as well as knock-down of Cfl1 in zebrafish result in morphological changes in the lower jaw and pharyngeal arches [96]. Our study showed that expression of Cfl1 was consistently higher in developing heads of the benthic morphs and its expression correlated negatively with Apaf1. Moreover, WISH analysis showed clear differences in the spatial expression pattern between the morphotypes, i.e. in the frontonasal region and in components of the lower jaw and pharyngeal arches (Fig. 7).

In zebrafish, an isoform of the transcription factor Foxq1 (Foxq1b) has been reported to be induced by activation of Ahr2 and its expression was limited to the jaw primordium [98]. Our results show higher expression of Foxq1 in developing heads of benthic morphs and strong negative expression correlation with Apaf1 (Figs. 5, 6a). We also identified differences in levels and spatial patterns of expression for Foxq1 between contrasting morphs in regions limited to ventral craniofacial structures, where it overlaps with the expression of Ahr2b, Apaf1, and Cfl1 (Fig. 7). Taken together, the transcriptional dynamics of Apaf1, Cfl1 and Foxq1 indicate higher AhR activity in benthic morphs and suggest a role for these three genes in divergent craniofacial morphogenesis in Arctic charr.

The other 5 AhR-regulated genes tested in this study, Cldn4, Jup, Mvp, Scin and Sult6b1, were found to be more highly expressed in benthic morphotypes (Fig. 5). Furthermore, most of these genes showed positive expression correlation with Ahr2b suggesting Ahr2b-dependent regulation. A junction plakoglobin gene, Jup, encodes a protein involved in cadherin-mediated intercellular communication [99]. In mammalian cells, Jup shows an interesting biphasic transcriptional response to AhR activity [100]. While high doses of AhR pathway agonists repressed Jup expression, a low dose could induce Jup expression [100]. Therefore, the increased expression of Jup in developing heads of benthic morphotypes might indicate moderate AhR pathway activity. Another interesting gene, Scin (Adseverin), was already shown in a previous transcriptome analysis to be among the top differentially expressed genes between two contrasting morphs [36]. Scin plays a role in rearranging the actin cytoskeleton, chondrocyte differentiation and in skeletal formation [101, 102]. In mice, Scin has been shown to be a direct target of the AhR pathway [103]. Additional studies are required to determine whether the genes found in this study have role in phenotypic variation of vertebrate craniofacial morphogenesis and whether and to what extent their transcription is regulated by the AhR pathway during skeletal development.

The differential transcriptional signatures of the selected AhR target genes clearly correlate with the differences in craniofacial development between charr morphotypes. However, further studies are required to differentiate between expression dynamics of the paralogues of each AhR target gene to identify the potential paralogue-specific response to the pathway. The quantified expression differences are most likely caused by genetic differences that affect regulators (TF’s, miRNAs, translation regulators or cis-elements) of the AhR pathway or parallel pathways. Can such differences in the developmental circuitry of craniofacial development play a role in the evolution of the charr morphs? Phenotypes and alleles of developmental genes can be influenced by genetic background [104, 105]. An interesting and highly relevant example is the fact that different inbred strains of mice display differential response to AhR-mediated effects on mandible development [106]. The morphs studied here all (except AC) exist in sympatry in Lake Thingvallavatn. Coalescent simulations involving PL and SB charr from the lake suggest a scenario of rapid divergence in allopatry followed by very slow (large populations and slow drift) divergence in sympatry and very low gene flow among the morphs [15]. Although the nature of the genetic separation (and possibly reproductive barriers) is unknown, our data on craniofacial development in hybrids of PL and SB charr suggested developmental incompatibilities resulting in smaller heads and narrowing of the mouth and pharyngeal tract in the hybrid offspring (Kapralova et al., unpublished data). This could reduce the fitness of hybrids (e.g. at the start of exogenous feeding), and lead to strong selection against hybrids. In this respect, it is tempting to speculate that the initial, rapid divergence processes postulated could involve changes in key developmental pathways of craniofacial development, possibly the AhR pathway, resulting in developmental incompatibility of hybrids, a situation that would select for further reproductive isolation during the sympatric phase.


In this study, we investigate the molecular mechanism underlying craniofacial divergence in the highly polymorphic Arctic charr of Lake Thingvallavatn (Iceland). We have found a network of co-expressed genes with higher expression in limnetic than benthic morphotypes during early stages of craniofacial development. Key members of the network and their genetic interactions are conserved across vertebrate species. Searching for predicted upstream TFs regulating the network, we identified a receptor of the AhR pathway, Ahr2b and found a set of AhR-regulated genes to be differentially expressed between the morphotypes. We also confirmed an overlapping expression pattern of Ahr2b and AhR target genes in craniofacial structures of Arctic charr embryos. Taken together, this study suggests the AhR pathway as a key modulator of transcriptional differences along the benthic–limnetic axis of craniofacial development in Arctic charr. Considering the conserved nature of the network, it is likely to have a much broader relevance in the early development of craniofacial structures in vertebrates.



aquaculture limnetic-like fish stock


quantification cycle


extracellular matrix


aryl hydrocarbon receptor pathway


large benthivorous charr


piscivorous charr


planktivorous charr


relative expression quantity


small benthivorous charr


transcription factor

τ s :

tau-somite relative age


  1. 1.

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

  2. 2.

    Helms JA, Cordero D, Tapadia MD. New insights into craniofacial morphogenesis. Development. 2005;132:851–61.

  3. 3.

    Yelick PC, Schilling TF. Molecular dissection of craniofacial development using zebrafish. Crit Rev Oral Biol Med. 2002;13:308–22.

  4. 4.

    Albertson RC, Kocher TD. Genetic and developmental basis of cichlid trophic diversity. Heredity (Edinb). 2006;97:211–21.

  5. 5.

    Ekblom R, Galindo J. Applications of next generation sequencing in molecular ecology of non-model organisms. Heredity (Edinb). 2011;107:1–15.

  6. 6.

    Wang S, Furmanek T, Kryvi H, Krossøy C, Totland GK, Grotmol S, Wargelius A. Transcriptome sequencing of Atlantic salmon (Salmo salar L.) notochord prior to development of the vertebrae provides clues to regulation of positional fate, chordoblast lineage and mineralisation. BMC Genom. 2014;15:141.

  7. 7.

    Braasch I, Peterson SM, Desvignes T, McCluskey BM, Batzel P, Postlethwait JH. A new model army: emerging fish models to study the genomics of vertebrate Evo-Devo. J Exp Zool B Mol Dev Evol. 2014;324:316–41.

  8. 8.

    Takechi M, Takeuchi M, Ota KG, Nishimura O, Mochii M, Itomi K, Adachi N, Takahashi M, Fujimoto S, Tarui H, Okabe M, Aizawa S, Kuratani S. Overview of the transcriptome profiles identified in hagfish, shark, and bichir: current issues arising from some nonmodel vertebrate taxa. J Exp Zool B Mol Dev Evol. 2011;316:526–46.

  9. 9.

    Skúlason S, Snorrason SS, Noakes DLG, Ferguson MM, Malmquist HJ. Segregation in spawning and early life history among polymorphic Arctic charr, Salvelinus alpinus, in Thingvallavatn, Iceland. J Fish Biol. 2006;35:225–32.

  10. 10.

    Snorrason SS, Skulason S, Jonsson B, Malquist HJ, Jonasson PM, Sandlund OT, Lindem T. Trophic specialization in Arctic charr Salvelinus alpinus (Pisces; Salmonidae): morphological divergence and ontogenetic niche shifts. Biol J Linn Soc. 1994;52:1–18.

  11. 11.

    Malmquist H, Snorrason S, Skúlason S, Jonsson B, Sandlund O, Jonasson P. Diet differentiation in polymorphic Arctic charr in Thingvallavatn, Iceland. J Anim Ecol. 1992;61:21–35.

  12. 12.

    Sandlund T, Gunnarsson K, Jónasson P, Jonsson B, Lindem T, Magnússon K, Malmquist H, Sigurjónsdóttir H, Skúlason S, Snorrason S. The Arctic Charr Salvelinus alpinus in Thingvallavatn. Oikos. 1992;64:305–51.

  13. 13.

    Jonsson B, Skúlason S, Snorrason SS, Sandlund OT, Malmquist HJ, Jónasson PM, Cydemo R, Lindem T. Life history variation of polymorphic Arctic Charr (Salvelinus alpinus) in Thingvallavatn, Iceland. Can J Fish Aquat Sci. 1988;45:1537–47.

  14. 14.

    Kapralova KH, Gudbrandsson J, Reynisdottir S, Santos CB, Baltanás VC, Maier VH, Snorrason SS, Palsson A. Differentiation at the MHCIIα and Cath2 loci in sympatric Salvelinus alpinus resource morphs in Lake Thingvallavatn. PLoS One. 2013;8:e69402.

  15. 15.

    Kapralova KH, Morrissey MB, Kristjánsson BK, Olafsdóttir GÁ, Snorrason SS, Ferguson MM. Evolution of adaptive diversity and genetic connectivity in Arctic charr (Salvelinus alpinus) in Iceland. Heredity (Edinb). 2011;106:472–87.

  16. 16.

    Skúlason S, Noakess D, Snorrason S. Ontogeny of trophic morphology in four sympatric morphs of arctic charr Salvelinus alpinus in Thingvallavatn, Iceland. Bioloqical J Linn Soc. 1989;38:281–301.

  17. 17.

    Parsons KJ, Skúlason S, Ferguson M. Morphological variation over ontogeny and environments in resource polymorphic arctic charr (Salvelinus alpinus). Evol Dev. 2010;12:246–57.

  18. 18.

    Küttner E, Parsons KJ, Robinson BW, Skúlason S, Danzmann RG, Ferguson MM. Effects of population, family, and diet on craniofacial morphology of Icelandic Arctic charr (Salvelinus alpinus). Biol J Linn Soc. 2013;108:702–14.

  19. 19.

    Ahi EP, Kapralova KH, Pálsson A, Maier VH, Gudbrandsson J, Snorrason SS, Jónsson ZO, Franzdóttir SR. Transcriptional dynamics of a conserved gene expression network associated with craniofacial divergence in Arctic charr. Evodevo. 2014;5:40.

  20. 20.

    Roberts RB, Hu Y, Albertson RC, Kocher TD. Craniofacial divergence and ongoing adaptation via the hedgehog pathway. Proc Natl Acad Sci USA. 2011;108:13194–9.

  21. 21.

    Parsons KJ, Trent Taylor A, Powder KE, Albertson RC. Wnt signalling underlies the evolution of new phenotypes and craniofacial variability in Lake Malawi cichlids. Nat Commun. 2014;5:3629.

  22. 22.

    Albertson RC, Yan Y-L, Titus TA, Pisano E, Vacchi M, Yelick PC, Detrich HW, Postlethwait JH. Molecular pedomorphism underlies craniofacial skeletal evolution in Antarctic notothenioid fishes. BMC Evol Biol. 2010;10:4.

  23. 23.

    Gunter HM, Koppermann C, Meyer A. Revisiting de Beer’s textbook example of heterochrony and jaw elongation in fish: calmodulin expression reflects heterochronic growth, and underlies morphological innovation in the jaws of belonoid fishes. Evodevo. 2014;5:8.

  24. 24.

    Powder KE, Milch K, Asselin G, Albertson RC. Constraint and diversification of developmental trajectories in cichlid facial morphologies. Evodevo. 2015;6:25.

  25. 25.

    Battle A, Khan Z, Wang SH, Mitrano A, Ford MJ, Pritchard JK, Gilad Y. Impact of regulatory variation from RNA to protein. Science (80). 2014;347:664–7.

  26. 26.

    Stuart JM, Segal E, Koller D, Kim SK. A gene-coexpression network for global discovery of conserved genetic modules. Science. 2003;302:249–55.

  27. 27.

    Schneider RF, Li Y, Meyer A, Gunter HM. Regulatory gene networks that shape the development of adaptive phenotypic plasticity in a cichlid fish. Mol Ecol. 2014;23:4511–26.

  28. 28.

    Filteau M, Pavey SA, St-Cyr J, Bernatchez L. Gene coexpression networks reveal key drivers of phenotypic divergence in lake whitefish. Mol Biol Evol. 2013;30:1384–96.

  29. 29.

    Levine M, Davidson EH. Gene regulatory networks for development. Proc Natl Acad Sci. 2005;102:4936–42.

  30. 30.

    Alon U. Network motifs: theory and experimental approaches. Nat Rev Genet. 2007;8:450–61.

  31. 31.

    Svavarsson E. Árangur í kynbótum á bleikju og næstu skref [reference in Icelandic]. In Fræðaþing landbúnaðarsins. 2007;4:121–5.

  32. 32.

    Gorodilov YN. Description of the early ontogeny of the Atlantic salmon, Salmo salar, with a novel system of interval (state) identification. Environ Biol Fishes. 1996;47:109–27.

  33. 33.

    Obayashi T, Kinoshita K. COXPRESdb: a database to compare gene coexpression in seven model animals. Nucleic Acids Res. 2011;39(Database issue):D1016–22.

  34. 34.

    Wang J, Duncan D, Shi Z, Zhang B. WEB-based GEne SeT AnaLysis Toolkit (WebGestalt): update 2013. Nucleic Acids Res. 2013;41:W77–83.

  35. 35.

    Kapralova KH, Jónsson ZO, Palsson A, Franzdóttir SR, Le Deuff S, Kristjanson BK, Snorrason SS. Bones in motion: ontogeny of craniofacial development in sympatric arctic charr morphs. Dev Dyn. 2015;244:1168–78.

  36. 36.

    Gudbrandsson J, Ahi E, Franzdottir S, Kapralova K, Kristjansson B, Steinhaeuser S, Maier V, Johannesson I, Snorrason S, Jonsson Z, Palsson A. The developmental transcriptome of contrasting Arctic charr (Salvelinus alpinus) morphs. F1000Research. 2015;4:136.

  37. 37.

    Moghadam HK, Ferguson MM, Danzmann RG. Whole genome duplication: challenges and considerations associated with sequence orthology assignment in Salmoninae. J Fish Biol. 2011;79:561–74.

  38. 38.

    Di Génova A, Aravena A, Zapata L, González M, Maass A, Iturra P. SalmonDB: a bioinformatics resource for Salmo salar and Oncorhynchus mykiss. Database (Oxford). 2011;2011:bar050.

  39. 39.

    Ramakers C, Ruijter JM, Deprez RHL, Moorman AF. Assumption-free analysis of quantitative real-time polymerase chain reaction (PCR) data. Neurosci Lett. 2003;339:62–6.

  40. 40.

    Ahi EP, Guðbrandsson J, Kapralova KH, Franzdóttir SR, Snorrason SS, Maier VH, Jónsson ZO. Validation of reference genes for expression studies during craniofacial development in Arctic Charr. PLoS One. 2013;8:e66389.

  41. 41.

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

  42. 42.

    Bergkvist A, Rusnakova V, Sindelka R, Garda JMA, Sjögreen B, Lindh D, Forootan A, Kubista M. Gene expression profiling–clusters of possibilities. Methods. 2010;50:323–35.

  43. 43.

    Team RDC. R: a language and environment for statistical computing. Austria: R Found Stat Comput Vienna; 2013.

  44. 44.

    Macqueen DJ, Johnston IA. Evolution of follistatin in teleosts revealed through phylogenetic, genomic and expression analyses. Dev Genes Evol. 2008;218:1–14.

  45. 45.

    Eroglu B, Wang G, Tu N, Sun X, Mivechi NF. Critical role of Brg1 member of the SWI/SNF chromatin remodeling complex during neurogenesis and neural crest induction in zebrafish. Dev Dyn. 2006;235:2722–35.

  46. 46.

    Bhattacherjee V, Horn KH, Singh S, Webb CL, Pisano MM, Greene RM. CBP/p300 and associated transcriptional co-activators exhibit distinct expression patterns during murine craniofacial and neural tube development. Int J Dev Biol. 2009;53:1097–104.

  47. 47.

    Sun F, Yang Q, Weng W, Zhang Y, Yu Y, Hong A, Ji Y, Pan Q. Chd4 and associated proteins function as corepressors of Sox9 expression during BMP-2-induced chondrogenesis. J Bone Miner Res. 2013;28:1950–61.

  48. 48.

    Kong Y, Grimaldi M, Curtin E, Dougherty M, Kaufman C, White RM, Zon LI, Liao EC. Neural crest development and craniofacial morphogenesis is coordinated by nitric oxide and histone acetylation. Chem Biol. 2014;21:488–501.

  49. 49.

    Rocha PP, Scholze M, Bleiss W, Schrewe H. Med12 is essential for early mouse development and for canonical Wnt and Wnt/PCP signaling. Development. 2010;137:2723–31.

  50. 50.

    Rau MJ, Fischer S, Neumann CJ. Zebrafish Trap230/Med12 is required as a coactivator for Sox9-dependent neural crest, cartilage and ear development. Dev Biol. 2006;296:83–93.

  51. 51.

    Ge W, Shi L, Zhou Y, Liu Y, Ma G-E, Jiang Y, Xu Y, Zhang X, Feng H. Inhibition of osteogenic differentiation of human adipose-derived stromal cells by retinoblastoma binding protein 2 repression of RUNX2-activated transcription. Stem Cells. 2011;29:1112–25.

  52. 52.

    Balow SA, Pierce LX, Zentner GE, Conrad PA, Davis S, Sabaawy HE, McDermott BM, Scacheri PC. Knockdown of fbxl10/kdm2bb rescues chd7 morphant phenotype in a zebrafish model of CHARGE syndrome. Dev Biol. 2013;382:57–69.

  53. 53.

    Bolli N, Payne EM, Rhodes J, Gjini E, Johnston AB, Guo F, Lee J-S, Stewart RA, Kanki JP, Chen AT, Zhou Y, Zon LI, Look AT. cpsf1 is required for definitive HSC survival in zebrafish. Blood. 2011;117:3996–4007.

  54. 54.

    Lines MA, Huang L, Schwartzentruber J, Douglas SL, Lynch DC, Beaulieu C, Guion-Almeida ML, Zechi-Ceide RM, Gener B, Gillessen-Kaesbach G, Nava C, Baujat G, Horn D, Kini U, Caliebe A, Alanay Y, Utine GE, Lev D, Kohlhase J, Grix AW, Lohmann DR, Hehr U, Böhm D, Majewski J, Bulman DE, Wieczorek D, Boycott KM. Haploinsufficiency of a spliceosomal GTPase encoded by EFTUD2 causes mandibulofacial dysostosis with microcephaly. Am J Hum Genet. 2012;90:369–77.

  55. 55.

    Quintana AM, Geiger EA, Achilly N, Rosenblatt DS, Maclean KN, Stabler SP, Artinger KB, Appel B, Shaikh TH. Hcfc1b, a zebrafish ortholog of HCFC1, regulates craniofacial development by modulating mmachc expression. Dev Biol. 2014;396:94–106.

  56. 56.

    Yamada-Okabe T, Imamura K, Kawaguchi N, Sakai H, Yamashita M, Matsumoto N. Functional characterization of the zebrafish WHSC1-related gene, a homolog of human NSD2. Biochem Biophys Res Commun. 2010;402:335–9.

  57. 57.

    Tuz K, Bachmann-Gagescu R, O’Day DR, Hua K, Isabella CR, Phelps IG, Stolarski AE, O’Roak BJ, Dempsey JC, Lourenco C, Alswaid A, Bönnemann CG, Medne L, Nampoothiri S, Stark Z, Leventer RJ, Topçu M, Cansu A, Jagadeesh S, Done S, Ishak GE, Glass IA, Shendure J, Neuhauss SCF, Haldeman-Englert CR, Doherty D, Ferland RJ. Mutations in CSPP1 cause primary cilia abnormalities and Joubert syndrome with or without Jeune asphyxiating thoracic dystrophy. Am J Hum Genet. 2014;94:62–72.

  58. 58.

    Wang B, Doan D, Roman Petersen Y, Alvarado E, Alvarado G, Bhandari A, Mohanty A, Mohanty S, Nissen RM. Wdr68 requires nuclear access for craniofacial development. PLoS One. 2013;8:e54363.

  59. 59.

    Caruana G, Bernstein A. Craniofacial dysmorphogenesis including cleft palate in mice with an insertional mutation in the discs large gene. Mol Cell Biol. 2001;21:1475–83.

  60. 60.

    Hall C, Flores MV, Murison G, Crosier K, Crosier P. An essential role for zebrafish Fgfrl1 during gill cartilage development. Mech Dev. 2006;123:925–40.

  61. 61.

    Iwase S, Lan F, Bayliss P, de la Torre-Ubieta L, Huarte M, Qi HH, Whetstine JR, Bonni A, Roberts TM, Shi Y. The X-linked mental retardation gene SMCX/JARID1C defines a family of histone H3 lysine 4 demethylases. Cell. 2007;128:1077–88.

  62. 62.

    Loomes KM, Stevens SA, O’Brien ML, Gonzalez DM, Ryan MJ, Segalov M, Dormans NJ, Mimoto MS, Gibson JD, Sewell W, Schaffer AA, Nah H-D, Rappaport EF, Pratt SC, Dunwoodie SL, Kusumi K. Dll3 and Notch1 genetic interactions model axial segmental and craniofacial malformations of human birth defects. Dev Dyn. 2007;236:2943–51.

  63. 63.

    Sutherland HG, Newton K, Brownstein DG, Holmes MC, Kress C, Semple CA, Bickmore WA. Disruption of Ledgf/Psip1 results in perinatal mortality and homeotic skeletal transformations. Mol Cell Biol. 2006;26:7201–10.

  64. 64.

    Wittschieben J, Shivji MKK, Lalani E, Jacobs MA, Marini F, Gearhart PJ, Rosewell I, Stamp G, Wood RD. Disruption of the developmentally regulated Rev3l gene causes embryonic lethality. Curr Biol. 2000;10:1217–20.

  65. 65.

    Goodale BC, La Du JK, Bisson WH, Janszen DB, Waters KM, Tanguay RL. AHR2 mutant reveals functional diversity of aryl hydrocarbon receptors in zebrafish. PLoS One. 2012;7:e29346.

  66. 66.

    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:4516–31.

  67. 67.

    Manousaki T, Hull PM, Kusche H, Machado-Schiaffino G, Franchini P, Harrod C, Elmer KR, Meyer A. 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:650–69.

  68. 68.

    Stock DW, Jackman WR, Trapani J. Developmental genetic mechanisms of evolutionary tooth loss in cypriniform fishes. Development. 2006;133:3127–37.

  69. 69.

    Shapiro MD, Marks ME, Peichel CL, Blackman BK, Nereng KS, Jónsson B, Schluter D, Kingsley DM. Genetic and developmental basis of evolutionary pelvic reduction in threespine sticklebacks. Nature. 2004;428:717–23.

  70. 70.

    Tanaka M, Hale LA, Amores A, Yan Y-L, Cresko WA, Suzuki T, Postlethwait JH. Developmental genetic basis for the evolution of pelvic fin loss in the pufferfish Takifugu rubripes. Dev Biol. 2005;281:227–39.

  71. 71.

    Albertson RC, Streelman JT, Kocher TD, Yelick PC. Integration and evolution of the cichlid mandible: the molecular basis of alternate feeding strategies. Proc Natl Acad Sci USA. 2005;102:16287–92.

  72. 72.

    Powder KE, Cousin H, McLinden GP, Craig Albertson R. A nonsynonymous mutation in the transcriptional regulator lbh is associated with cichlid craniofacial adaptation and neural crest cell development. Mol Biol Evol. 2014;31:3113–24.

  73. 73.

    Fraser GJ, Hulsey CD, Bloomquist RF, Uyesugi K, Manley NR, Streelman JT. An ancient gene network is co-opted for teeth on old and new jaws. PLoS Biol. 2009;7:e31.

  74. 74.

    Knight RD, Javidan Y, Zhang T, Nelson S, Schilling TF. AP2-dependent signals from the ectoderm regulate craniofacial development in the zebrafish embryo. Development. 2005;132:3127–38.

  75. 75.

    Barouki R, Coumoul X, Fernandez-Salguero PM. The aryl hydrocarbon receptor, more than a xenobiotic-interacting protein. FEBS Lett. 2007;581:3608–15.

  76. 76.

    Nguyen LP, Bradfield CA. The search for endogenous activators of the aryl hydrocarbon receptor. Chem Res Toxicol. 2008;21:102–16.

  77. 77.

    Casado FL, Singh KP, Gasiewicz TA. The aryl hydrocarbon receptor: regulation of hematopoiesis and involvement in the progression of blood diseases. Blood Cells Mol Dis. 2010;44:199–206.

  78. 78.

    Sartor MA, Schnekenburger M, Marlowe JL, Reichard JF, Wang Y, Fan Y, Ma C, Karyala S, Halbleib D, Liu X, Medvedovic M, Puga A. Genomewide analysis of aryl hydrocarbon receptor binding targets reveals an extensive array of gene clusters that control morphogenetic and developmental programs. Environ Health Perspect. 2009;117:1139–46.

  79. 79.

    Fujii-Kuriyama Y, Kawajiri K. Molecular mechanisms of the physiological functions of the aryl hydrocarbon (dioxin) receptor, a multifunctional regulator that senses and responds to environmental stimuli. Proc Japan Acad Ser B. 2010;86:40–53.

  80. 80.

    Anderson G, Beischlag TV, Vinciguerra M, Mazzoccoli G. The circadian clock circuitry and the AHR signaling pathway in physiology and pathology. Biochem Pharmacol. 2013;85:1405–16.

  81. 81.

    Moura-Alves P, Faé K, Houthuys E, Dorhoi A, Kreuchwig A, Furkert J, Barison N, Diehl A, Munder A, Constant P, Skrahina T, Guhlich-Bornhof U, Klemm M, Koehler A-B, Bandermann S, Goosmann C, Mollenkopf H-J, Hurwitz R, Brinkmann V, Fillatreau S, Daffe M, Tümmler B, Kolbe M, Oschkinat H, Krause G, Kaufmann SHE. AhR sensing of bacterial pigments regulates antibacterial defence. Nature. 2014;512:387–92.

  82. 82.

    Liu Q, Spitsbergen JM, Cariou R, Huang C-Y, Jiang N, Goetz G, Hutz RJ, Tonellato PJ, Carvan MJ. Histopathologic alterations associated with global gene expression due to chronic dietary TCDD exposure in juvenile zebrafish. PLoS One. 2014;9:e100910.

  83. 83.

    Kung T, Murphy KA, White LA. The aryl hydrocarbon receptor (AhR) pathway as a regulatory pathway for cell adhesion and matrix metabolism. Biochem Pharmacol. 2009;77:536–46.

  84. 84.

    Kel A, Reymann S, Matys V, Nettesheim P, Wingender E, Borlak J. A novel computational approach for the prediction of networked transcription factors of aryl hydrocarbon-receptor-regulated genes. Mol Pharmacol. 2004;66:1557–72.

  85. 85.

    Hansson MC, Hahn ME. Functional properties of the four Atlantic salmon (Salmo salar) aryl hydrocarbon receptor type 2 (AHR2) isoforms. Aquat Toxicol. 2008;86:121–30.

  86. 86.

    Prasch AL, Teraoka H, Carney SA, Dong W, Hiraga T, Stegeman JJ, Heideman W, Peterson RE. Aryl hydrocarbon receptor 2 mediates 2,3,7,8-tetrachlorodibenzo-p-dioxin developmental toxicity in zebrafish. Toxicol Sci. 2003;76:138–50.

  87. 87.

    Xiong KM, Peterson RE, Heideman W. Aryl hydrocarbon receptor-mediated down-regulation of sox9b causes jaw malformation in zebrafish embryos. Mol Pharmacol. 2008;74:1544–53.

  88. 88.

    Huang Z, Xu H, Sandell L. Negative regulation of chondrocyte differentiation by transcription factor AP-2alpha. J Bone Miner Res. 2004;19:245–55.

  89. 89.

    Iwata H, Stegeman JJ. In situ RT-PCR detection of CYP1A mRNA in pharyngeal epithelium and chondroid cells from chemically untreated fish: involvement in vertebrate craniofacial skeletal development? Biochem Biophys Res Commun. 2000;271:130–7.

  90. 90.

    Carney SA, Peterson RE, Heideman W. Hydrocarbon receptor/aryl hydrocarbon receptor nuclear translocator pathway causes developmental toxicity through a CYP1A-independent mechanism in zebrafish. Mol Pharmacol. 2004;66:512–21.

  91. 91.

    Cecconi F, Alvarez-bolado G, Meyer BI, Roth KA, Gruss P, Gene A. Apaf1 (CED-4 homolog) regulates programmed cell death in mammalian development generation of a null allele of the murine. Cell. 1998;94:727–37.

  92. 92.

    Long AB, Kaiser WJ, Mocarski ES, Caspary T. Apaf1 apoptotic function critically limits Sonic hedgehog signaling during craniofacial development. Cell Death Differ. 2013;20:1510–20.

  93. 93.

    Hu Y, Albertson RC. Hedgehog signaling mediates adaptive variation in a dynamic functional system in the cichlid feeding apparatus. Proc Natl Acad Sci. 2014;111:8530–4.

  94. 94.

    Marlowe JL, Fan Y, Chang X, Peng L, Knudsen ES, Xia Y, Puga A. The aryl hydrocarbon receptor binds to E2F1 and inhibits E2F1-induced apoptosis. Mol Biol Cell. 2008;19:3263–71.

  95. 95.

    Marlowe JL, Knudsen ES, Schwemberger S, Puga A. The aryl hydrocarbon receptor displaces p300 from E2F-dependent promoters and represses S phase-specific gene expression. J Biol Chem. 2004;279:29013–22.

  96. 96.

    Ashworth S, Teng B, Kaufeld J, Miller E, Tossidou I, Englert C, Bollig F, Staggs L, Roberts ISD, Park J-K, Haller H, Schiffer M. Cofilin-1 inactivation leads to proteinuria–studies in zebrafish, mice and humans. PLoS One. 2010;5:e12626.

  97. 97.

    Huang L, Wang C, Zhang Y, Li J, Zhong Y, Zhou Y, Chen Y, Zuo Z. Benzo[a]pyrene exposure influences the cardiac development and the expression of cardiovascular relative genes in zebrafish (Danio rerio) embryos. Chemosphere. 2012;87:369–75.

  98. 98.

    Planchart A, Mattingly CJ. 2,3,7,8-Tetrachlorodibenzo-p-dioxin upregulates FoxQ1b in zebrafish Jaw primordium. Chem Res Toxicol. 2010;23:480–7.

  99. 99.

    Cowin P, Kapprell H-P, Franke WW, Tamkun J, Hynes RO. Plakoglobin: A protein common to different kinds of intercellular adhering junctions. Cell. 1986;46:1063–73.

  100. 100.

    Procházková J, Kozubík A, Machala M, Vondráček J. Differential effects of indirubin and 2,3,7,8-tetrachlorodibenzo-p-dioxin on the aryl hydrocarbon receptor (AhR) signalling in liver progenitor cells. Toxicology. 2011;279:146–54.

  101. 101.

    Vieira FA, Thorne MAS, Stueber K, Darias M, Reinhardt R, Clark MS, Gisbert E, Power DM. Comparative analysis of a teleost skeleton transcriptome provides insight into its regulation. Gen Comp Endocrinol. 2013;191:45–58.

  102. 102.

    Nurminsky D, Magee C, Faverman L, Nurminskaya M. Regulation of chondrocyte differentiation by actin-severing protein adseverin. Dev Biol. 2007;302:427–37.

  103. 103.

    Svensson C, Silverstone AE, Lai Z-W, Lundberg K. Dioxin-induced adseverin expression in the mouse thymus is strictly regulated and dependent on the aryl hydrocarbon receptor. Biochem Biophys Res Commun. 2002;291:1194–200.

  104. 104.

    Chari S, Dworkin I. The conditional nature of genetic interactions: the consequences of wild-type backgrounds on mutational interactions in a genome-wide modifier screen. PLoS Genet. 2013;9:e1003661.

  105. 105.

    Chandler CH, Chari S, Dworkin I. Does your gene need a background check? How genetic background impacts the analysis of mutations, genes, and evolution. Trends Genet. 2013;29:358–66.

  106. 106.

    Keller JM, Zelditch ML, Huet YM, Leamy LJ. Genetic differences in sensitivity to alterations of mandible structure caused by the teratogen 2,3,7,8-tetrachlorodibenzo-p-dioxin. Toxicol Pathol. 2008;36:1006–13.

Download references

Authors’ contributions

EPA designed and performed most experiments, conducted data analysis and drafted the manuscript. SSSt assisted with literature studies, selection of gene candidates, design of primers and performed some qPCR experiments. AP assisted in statistical data analysis. VHM set up in situ hybridization experiments. AP, SRF, SSSn, VHM and ZOJ designed the study, acquired funding, interpreted the data and co-wrote the manuscript. All authors read and approved the final manuscript.


This project was supported by The Icelandic Centre for Research (RANNIS/IRF, Grant 100204), The University of Iceland Research Fund and the Eimskip University Fund. The authors thank Bjarni K Kristjánsson, Einar Svavarsson, Kalina H. Kapralova and Soizic Le Deuff for their assistance with sampling the parents, setting up crosses, maintenance and sampling of the embryos and experimental advice and discussion. The authors also thank Jóhannes Gudbrandsson for providing R-code for cluster analysis of gene expression data, statistical advice and helpful discussions. Finally, we thank Eggert Pétursson for permitting us to use his drawings of Arctic charr heads.

Compliance with ethical guidelines

Competing interests The authors declare that they have no competing interests.

Author information

Correspondence to Valerie H. Maier.

Additional files


Additional file 1: Table S1. AhR-regulated candidate genes and available literature on their response to AhR activation, expression during craniofacial development and involvement in skeletal formation.


Additional file 2: Table S2. Primers designed for cDNA cloning and qPCR study of Arctic charr.


Additional file 3: Table S3. Chd4-Cspp1-Dlg1 co-expressed genes in humans and mice retrieved from COXPRESdb. More than 200 genes were identified to have co-expression with Chd4, Cspp1 and Dlg1 after setting the mutual rank (MR) to the top-ranked 2000 genes and the reliability score of 1 to 3.


Additional file 4: Table S4. Transcription factor enrichment analysis using WEB-based GEne SeT AnaLysis Toolkit (WebGestalt v2). The list of more than 200 genes co-expressed with Chd4, Cspp1 and Dlg1 in human and mice are used as input. (Cut-off criteria: p < 0.001 adjusted with hypergeometric statistical analysis; Bonferroni multiple test adjustment).


Additional file 5: Figure S1. Analyses of correlations between the expression of TFs (E2f1, Yy1, Nrf1, Sp1, Maz and Myc) and expression of the candidate target genes in the study classified in three groups. The Up group indicates genes with significantly higher expression in the benthic charr morphs, Down includes genes with significantly lower expression in benthic morphs and Rest constitutes genes without significant expression difference between benthic and limnetic charr. The boxplots summarize Kendall correlation coefficients by gene groups. Non parametric Kruskal–Wallis ANOVA indicated differences among gene groups for all TFs except Myc (P = 0.015 for Sp1, P < 0.005 for the other four genes).

Rights and permissions

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

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark


  • Arctic Charr
  • Aryl Hydrocarbon Receptor
  • Pharyngeal Arch
  • Craniofacial Development
  • Aryl Hydrocarbon Receptor Pathway