The developmental transcriptome of contrasting Arctic charr ( Salvelinus alpinus) morphs

Species and populations with parallel evolution of specific traits can help illuminate how predictable adaptations and divergence are at the molecular and developmental level. Following the last glacial period, dwarfism and specialized bottom feeding morphology evolved rapidly in several landlocked Arctic charr Salvelinus alpinus populations in Iceland. To study the genetic divergence between small benthic morphs and limnetic morphs, we conducted RNA-sequencing charr embryos at four stages in early development. We studied two stocks with contrasting morphologies: the small benthic (SB) charr from Lake Thingvallavatn and Holar aquaculture (AC) charr. The data reveal significant differences in expression of several biological pathways during charr development. There was also an expression difference between SB- and AC-charr in genes involved in energy metabolism and blood coagulation genes. We confirmed differing expression of five genes in whole embryos with qPCR, including lysozyme and natterin-like which was previously identified as a fish-toxin of a lectin family that may be a putative immunopeptide. We also verified differential expression of 7 genes in the developing head that associated consistently with benthic v.s.limnetic morphology (studied in 4 morphs). Comparison of single nucleotide polymorphism (SNP) frequencies reveals extensive genetic differentiation between the SB and AC-charr (~1300 with more than 50% frequency difference). Curiously, three derived alleles in the otherwise conserved 12s and 16s mitochondrial ribosomal RNA genes are found in benthic charr. The data implicate multiple genes and molecular pathways in divergence of small benthic charr and/or the response of aquaculture charr to domestication. Functional, genetic and population genetic studies on more freshwater and anadromous populations are needed to confirm the specific loci and mutations relating to specific ecological traits in Arctic charr.

Historical contingencies and chance shape organisms during evolution 1,2 , but convergence in phenotype and molecular systems indicates that evolution is to some extent predictable 3,4 . Identification of genes and variants that influence evolved differences is not a trivial task 5 . Ideal systems to study the role of chance and necessity in ecological evolution would be related species or populations with readily observable phenotypic variation, living in a tractable ecological setting, and showing parallel evolution of specific traits within/among species/populations. Examples of such species complexes are provided finches of the Galapagos islands 6 , while cichlids of the African great lakes also provide an exciting mulit-species system in the same respect 7 . The threespine stickleback has also emerged as a model "single species" system 8 . The amount of diversity in the feeding specializations of fish provide great opportunities for studying adaptation and divergence at the developmental and genetic level.
One approach to identify pathways related to function or morphological differences between species, populations or ecomorphs is to study gene expression during development 9,10 . For example a microarray study of liver samples from anadromous and resident populations of brown trout (Salmo trutta), revealed that gene expression in juveniles was more influenced by life history than relatedness 11 . Furthermore, Filteau et al. (2013) 12 found a set of coexpressed genes differenting two whitefish morphotypes, implicating Bone morphogenesis protein (BMP) signaling in the development of ecological differences in trophic morphology. Thus we were quite keen to apply RNA-sequencing to analyze ecomorphs in our study system, Arctic charr. Two previous studies have used RNA-seq to study salinity tolerance in adult Arctic charr, and found links between gene expression and quantitative trait loci 13,14 .
Some northern freshwater fish species exhibit frequent parallelism in trophic structures and life history and in several cases are found as distinct resource morphs 8,[15][16][17][18][19] . One of these species, Arctic charr (Salvelinus alpinus), is well suited for studying the developmental underpinnings of trophic divergence and parallel evolution. Local adaptation has been extensively studied in the salmonid family, to which Arctic charr belongs 20 . The family is estimated to be between 88-103 million years old 21,22 . A whole genome duplication event occurred before the radiation of the salmonid family 21-24 which has provided time for divergence of ohnologous genes (paralogous genes originated by whole genome duplication event). Furthermore, recent estimates from the rainbow trout (Oncorhynchus mykiss) genome suggest that ohnologous genes are lost at a rate of about 170 genes per million years, and that on the order of 4500 were retained in rainbow trout 22 . De novo assembly of genomes and transcriptomes is complicated if many paralogs are present, which is the case in Arctic charr -see 13,14. In this study we opted for mapping the reads (36 bp) to a related reference genome/transcriptome 25 , instead of de novo assembly.
Molecular studies of the highly polymorphic Arctic charr Following the end of the last glacial period, about 10.000 years ago, Arctic charr colonized northern freshwater systems 26 . It is found as anadromous or lake/stream residents and exhibits high level of within species polymorphism 17,26 . Charr is also known to harbour substantial phenotypic plasticity, which may promote or reduce divergence 15,27 . Resource polymorphism in charr correlates with ecological attributes [28][29][30] . For instance small charr with benthic morphology, are found in multiple lavaspring and pond habitats in Iceland 31 , and a comparative study of Icelandic lakes 30 found that lakes with greater limnetic habitat, lower nutrients levels, and greater potential for zooplankton consumption appeared to promote resource polymorphism. Some of the larger lakes contain two or more distinct morphs, typically limnetic and benthic forms. Multiple lines of evidence show that these differences stem both from environmental and genetic causes 32-36 . The best studied example of sympatric charr are the four morphs in Lake Thingvallavatn 37 ; two have a benthic morphotype, a large benthivorous (LB-charr) and a small benthivorous (SB-charr), and two morphs are limnetic, a large piscivorous morph (PI-charr) and small planktivorous morph (PL-charr) 38 . Both PL and PI-charr operate in open water and feed on free-swimming prey, PL on planktonic crustaceans and PI on small fish. The PL, LB and SB-charr are presented in Figure 1.
Several population genetics studies, using allozymes or mtDNA revealed no differences among charr morphs in Lake Thingvallavatn 39-41 while other studies using microsatellite markers and nuclear genes, found significant 42-44 genetic differences among morphs in the lake 45 . Importantly Kapralova et al. (2011) 44 concluded that small benthic morphs have evolved repeatedly in Iceland and that gene flow has been reduced between the PL and SB morphs in Lake Thingvallavatn since its formation approximately 10,000 years ago 46 . We also discovered genetic separation in immunological genes (MHCIIα and cath2) between morphs in Iceland and within the lake 45 , consistent with ecologically driven evolution of immune functions. Recently qPCR analyses showed that expression of mTOR pathway components in skeletal muscle correlates with the SB-charr form in Iceland 47 , but it is unknown whether there is genetic differentiation in those genes or upstream

Amendments from Version 1
The major changes to the manuscript involve rewriting of the introduction, to highlight the differences between the overall objective of our research program and the objectives of this study. We are interested in studying the genetics of parallel evolution, but this study focuses on revealing differences in expression and genes separating sympatric benthic and limnetic morphs. We added clarifications on several aspects of the work and analyses, for instance providing workflow and sample overview in new Figure 2, adding morphs and descriptions to Figure 1, explaining the sampling and SNP filtering.
The reviewers pointed out a mistake in our interpretation of the fate of paralogous genes in salmonids, based on the Rainbow trout genome, which we met by fixing the manuscript and changing the interpretations. The reviewers questioned the choice of transcripts studied with qPCR. We rewrote this section, and added a table that summarizes the qPCR data, and highlights the fact that the data can be used to find candidate genes or paralog groups with expression differences between Arctic charr morphs. The reviewers also highlighted the low efficiency in the qPCR reactions on the natterin-like paralogs. We tried to accommodate those weaknesses by redoing figures, with the appropriate correction factors.
Multiple other aspects of the text where rewritten in accordance with the suggestions of the reviewers, which have in our opinion greatly improved the manuscript.

REVISED
regulators. Because individual genes have distinct histories 48,49 , genome wide methods are needed to identify genes and mutation that associate with divergence. Icelandic aquaculture charr (AC) was founded with fish from the north of Iceland, and has been bred at Holar University College since 1990 50 . The Holar AC-charr has responded to artificial selection in growth and performance characteristics, and is now the dominant charr breed in aquaculture in Iceland. While clearly a derived form, it has retained general limnetic craniofacial morphotype ( Figure 1).
In this study we compare SB-charr from Lake Thingvallavatn and AC-charr because i) SB charr represents an extensively studied and derived form of charr, that has been separated from Anadromous fish for approx. 10,000 years, ii) of the availability of abundant AC material and iii) we wanted an extreme contrast, because of budget reasons we could only sequence 8 samples at the time. The AC-charr is included here as a limnetic reference population, in part because we were unable to catch spawning anadromous charr, the ideal outgroup. But by focusing the follow up work on sympatric benthic and limnetic morphs of Lake Thingvallavatn, we can test and verify a subset of the signals found here. The contrast of SB and AC is justified as the data and studies (51-53) building on this data illustrate (see discussion).
The overall objectives of our research program are to investigate the genetics and developmental underpinnings of charr divergence and benthic parallelism. As a step towards this we compare the developmental transcriptome of SB charr and AC charr. The aims of this study are threefold. First, to find genes and pathways related to the development of phenotypic differences between small benthic charr from Lake Thingvallavatn and Icelandic aquaculture charr conforming to a limnetic morphotype. Second, to screen for signals of genetic differentiation between these two charr types. Third, we set out to verify a subset of the expression and genetic signals, in these morphs and two more (benthic and limnetic) morphs from Lake Thingvallavatn. We conduct RNA-sequencing of developing offspring of these two contrasting Arctic charr morphs, reared in common lab environment to minimize the effects of environmentally Adult individuals of the four morphs studied here, from above; the Holar aquaculture charr, the small benthic charr; the planktivorous charr; and the large benthic charr. The latter three all come from Lake Thingvallavatn and were sexually ripe. The morphs differ in size at maturation, body and head shape -mainly lower jaw and length of maxilla and colour pattern in the wild.
induced phenotypic plasticity on developmental phenotypes and gene expression. The data reveal genetic changes in nuclear and mitochondrial genes and differential expression of genes that may affect craniofacial and phenotypic traits which separate benthic and limnetic morphotypes in charr.

Methods
Sampling, rearing and developmental series Overview of the experimental design, RNA sequencing, analyses and follow work is outlined in Figure 2. We set up crosses and reared embryos in the laboratory as described in 51. Embryos from four charr morphs were studied: an aquaculture charr (AC-charr) from the Holar breeding program 50 and three natural morphs from Lake Thingvallavatn; SB, LB and PL-charr 54 . Samples of the first two, AC and SB-charr, with contrasting adult size and morphology (Figure 1), were collected in 2009 and material sent for RNA sequencing. The latter two were sampled in 2010 and were used for qPCR and SNP studies of selected genes. Briefly, in September 2009 we got material from spawning AC-charr from the Holar breeding program, from single parent crosses 50 and spawning SB-charr collected via gill netting in Olafsdrattur in Lake Thingvallavatn. Similarly, in the 2010 spawning season SB-, LB-and PL-charr were collected from Lake Thingvallavatn. For each parent group, eggs from several females (3-10) were pooled and fertilized using milt from several males (3-5) from the same group. Embryos were reared at ~ 5°C under constant water flow and in complete darkness at the Holar University College experimental facilities in Verid, Saudárkrókur. The water temperature was recorded twice daily and the average was used to estimate the relative age of the embryos using tausomite units (τs) 55 . Embryos and juveniles were sampled at designated time points, placed in RNAlater (Ambion) and frozen at −20°C. Post hatching juveniles were reared at the same temperature on standard Aquaculture food. For the investigation of different tissues of adult aquaculture charr (AC) from Hólar (fish size 20-25 cm) were used. Six randomly selected individuals were killed (by cutting through spinal cord) and dissected, and samples were taken from the skin, heart, liver, gills, spleen, intestine and kidney of each fish. The samples were placed in RNAlater (Ambion) and stored at −20°C. We used DNA for population genetic analyses from our previous study 45 , eight individuals from each of the three types, PL, LB and SB-charr.
Fishing in Lake Thingvallavatn was done with permissions obtained both from the owner of the land in Mjóanes and from the Thingvellir National Park commission. Ethics committee approval is not needed for regular or scientific fishing in Iceland (The  Icelandic law  To estimate expression, reads were aligned with RSEM version 1.1.18 with default parameters. RSEM distributes reads that map to multiple locations to the most likely contig, using expectation maximization 58 . The read counts for contigs with the same annotation were pooled because some genes were represented by more than one contig, and due to whole genome duplication almost the half of salmonid genes exist as ohnologs 22,24 . Thus the expression tests are done on gene or paralog group level, instead of the contig level. We acknowledge that paralogous genes are not always expressed similarly, but feel its necessary to do this pooling because of the nature of the data. In the remainder of the paper, we will refer to gene or paralog group (the number of underlying contigs is indicated in relevant tables). This brought the number of genes considered down to 16851. Lastly, paralog groups with fewer than 800 mapped reads in the entire dataset were excluded from the analyses, yielding a total of 10496.
A generalized linear model (GLM) with morph and developmental time as explanatory variables was used to find genes with different expression levels between the two charr morphotypes (groups) using the edgeR-package in R 59 .

Y = Morph + Time + Error
To obtain further insight into the expression profiles of differently expressed genes, we performed clustering analyses on log-transformed cpm-values (counts per million; cpm-function in edgeR). The values for each gene were scaled by mean and standard deviation, and the euclidean distance used for the hclust-function in R 60 with the default settings. We used the hypergeometric-test in goseq 61 to test for gene ontology enrichment. Since we pooled the read-count from different contigs we could unfortunately not take gene length into account in those tests.

Tests of differential expression with qPCR
We previously identified suitable reference genes to study Arctic charr development 51 . Here we examined the expression of several genes in whole charr embryos, embryonic heads and adult tissues. Primers were designed using the Primer3 tool 62 and checked for self-annealing and heterodimers according to the MIQE guidelines 63 (S1 Table). Primers for genes with several paralogs were designed for regions conserved among paralogs, except for natterin-like, where primers were designed to match regions differing in sequence between paralogs. Relative expression was calculated using the 2 −ΔΔCt method 64 . For the calculation of relative expression of genes in whole embryos, the geometric mean expression of three reference genes, β-Actin (Actb), elongation factor 1α and Ubiquitin-conjugating enzyme E2 L3, was used for normalization. For visual comparisons among samples, the normalized expression was presented as relative to the expression in AC at 161 τ s (calibration sample). For the embryonic head samples Eukaryotic Translation Initiation Factor 5A (If5a1) and Actb were used as reference genes and a biological replicate of AC at 178 (τs ) as the calibrator sample, see 51,52. Standard errors of relative expression were calculated from the standard errors (SE) of the ΔC T -values with the formula 2 −(ΔΔCt+SE) = minimum fold expression and 2 −(ΔΔCt−SE) = maximum fold expression. The statistical analysis was performed using the ΔC T -values with a two-way ANOVA with GLM function in R.

Y = Morph + Time + M x T + Error
Normal distribution of residuals was confirmed for all data. For the study of expression in the embryonic head we followed a significant morph effect in the ANOVA with Tukey's post-hoc honest significant difference test, on relative expression ratios (ΔC T s). Three genes had lower efficiency (as low as 1.72). We acknowledge that the data on those genes may be weak.

Polymorphisms in the Arctic charr transcriptome
For analysis of genetic variation we mapped the reads to the salmon contigs, this time using the Burrow-Wheeler Aligner (BWA) 65 with a seed length of 25, allowing two mismatches. We re-mapped the reads, since BWA allows short indels (RSEM does not) but disregarding them leads to many false SNPs close to indels. To extract candidate polymorphic sites from the Arctic charr transcriptome we ran VarScan2 66 with minimum coverage of 50 reads and minimum minor allele frequency of 0.1 on reads mapped to each S. salar contig for all of the 8 timepoints and morph combinations. This was done separately for reads that mapped uniquely to one contig only (UNI) and reads that mapped to two or more contigs (REP). These SNP-candidates were further processed in R 60 , following established principles for variant calling 67 . SNP-candidates at 90% frequency or higher in all samples were disregarded, as they reflect differences between Arctic charr and S. salar and are not the focus of this study. SNP-candidates with poor coverage in specific samples -i.e. coverage of five or fewer reads in three or four samples of each morph -were removed. As the SNP analysis was done on individual contigs, differences among paralogs appear in the data.
To address this we use the fact that each sample is a pool of few individuals, thus true SNPs are unlikely to have the same frequency in all samples. However, variants reflecting differences between paralogs will have similar frequency all samples (assuming steady difference in their expression in all samples). We evaluated differences between samples with Fisher exact tests, and only SNPs significantly different between samples with a p < 0.05 (with no multiple testing correction) were retained. To compare morphs, read numbers were summed over the four samples from each morph. A conservative approach was taken by focusing on SNP-candidates that showed the largest differences in frequency between morphs (delta), without adjusting for multiple testing (Fisher exact test, p > 5%). SNP-candidates with the highest frequency difference (delta > 95%) were manually processed and redundant candidates removed. A similar approach was used to mine for polymorphisms in Arctic charr mtDNA (NC_000861), using S. salar mtDNA as the outgroup (NC_001960.1).
We wrote a python script to predict the impact of SNPs within the mRNA sequences. Polymorphisms were categorized according to their location (3'UTR, coding, 5'UTR), and those within the coding region into synonymous or non-synonymous.

Verification of candidate SNPs
We chose 12 candidate SNPs for verification (see below). As the AC-charr is not a random breeding population, and because our interest is on differences between wild morphs, we took random samples of spawning SB, LB and PL-charr from Lake Thingvallavatn (8 per morph) from our earlier study 45 . Using the same PCR and DNA sequencing approach we genotyped 12 candidate SNPs (S2 Table).

RNA sequencing characteristics
Each sample yielded good quality data, with sequencing depth from 49 to 58 million (average: 55 million) reads. To quantify the expression levels, the reads were aligned to a salmon EST-assembly 57 . Around 20% of the reads mapped uniquely to the EST data (S3 Table). A further 30% mapped to two or more contigs, probably representing paralogous genes, recent duplications or repeat-like elements within transcribed regions. A substantial fraction of the RNA-sequencing reads did not map to the contigs from S. salar. Analyses of those reads require an Arctic charr genome sequence or transcriptome assembly from longer and paired end reads, currently underway in our laboratory.

Differential expression during Arctic charr development
We detected considerable changes in the transcriptome during Arctic charr development ( Figure 3a). The expression of 1603 and 2459 paralog groups differed significantly between developmental timepoints at the 1% and 5% levels of false discovery rate (FDR), respectively (Dataset 1). The difference was most pronounced between prehatching (timepoints: 141, 163, 200 τs) and post hatching embryos (timepoint 433 τs), as more than 70% of the paralog groups with FDR below 1% had higher expression in the latter ( Figure 3b). Gene Ontology analyses reveal six enriched GO categories (below 10%FDR). The most drastic changes were seen in processes related to glycolysis (GO:0006096, FDR = 0.0009), where the expression of 19 out of 25 paralog groups changed during this developmental period. The other five classes that were differentially expressed during charr development are: ion transport (GO:0006811, FDR = 0.027), blood coagulation (GO:0007596, FDR = 0.03), DNA repair (GO:0006281, FDR = 0.08) and two immune related categories (GO:0019882, FDR = 0.08, GO:0006955, FDR = 0.09). Those results probably reflect developmental changes and/or differences in the environment of embryos before and after hatching.

Differential expression between Arctic charr morphs
We were especially interested in genes showing expression differences between the two morphs as they might implicate pathways involved in the ecological divergence among charr populations.
In the data 296 paralog groups were differentially expressed (FDR < 5%) between the morphs (141 higher in SB and 152 higher in AC-charr, Dataset 1). Among genes with higher expression in SB-charr two biological GO categories were enriched: blood coagulation (GO:0007596, p = 0.001) and proteolysis (GO:0006508, p = 0.002). Recall, expression of blood coagulation factors also differed between developmental stages (see above). In AC-charr, genes in three categories: respiratory electron transport chain (GO:0022904, p = 0.0006), ATP synthesis coupled electron transport (GO:0042773, p = 0.002) and neurotransmitter transport (GO:0006836, p = 0.009) have higher expression. The first two GO categories both relate to energy generation in mitochondria and could reflect higher expression of genes with mitochondrial functions in AC-charr. At more stringent FDR (1%), 31 paralog groups, with diverse functional annotations, were higher expressed in SB and 40 genes higher in AC-charr ( Figure 3b, Table 1 and Table 2). The higher expressed ESTs were clustered into 4 groups for each morph, reflecting in some cases functional similarity. For instance SB cluster 3 has three immune related paralog groups: Complement factor D (9), H-2 class I histocompatibility antigen L-D alpha chain (2) and Sushi domain-containing protein 2 (4) ( Table 1). Note, however, that immune genes were not significantly enriched in the GO comparison of morphs. The results suggest genes with mitochondrial function, blood coagulation and other functions are differentially expressed between the morphs, but as few samples were sequenced, qPCR verification was needed.
Validation of gene expression differences in whole embryos and paralog specific expression of natterin genes For validation we opted for qPCR analyses of 9 genes/paralog groups in whole embryos and 8 in embryonic heads, which showed differential expression between AC and SB-charr, with statistical support ranging from <1% to about 10% FDR. We studied paralog groups with less FDR support, in part to be able to cast a wider net (see below). Of the nine paralog groups studied in whole embryos, five were confirmed to be differentially expressed between AC     Table and Dataset 2). Three genes, Nattl, Alkaline phosphatase (Alp) and Lysozyme C II (Lyz2), had significantly higher expression in SB. The other two, Keratin-associated protein 4-3 (Krtap4-3) and Poly polymerase 6 (Parp6) had higher expression in AC embryos ( Figure 4, S4 Table). No morph and time interaction was detected for any of the genes.
As some genes are represented by different contigs or even paralogs, we set out to disentangle the expression of one paralog group (Nattl) in detail. We measured the expression of three natterin paralogs (nattl1, nattl2 and nattl3), by designing qPCR primers that matched divergent regions. These genes caught our interest because the only prior work implicated Natterin as a toxin produced by a tropical fish 70,71 . We studied nattl expression in several developmental stages in AC-, SB-and PL-charr as well as in selected tissues of adult AC-charr. The expression level of the three paralogs differed between morphs and timepoints ( Figure 5 and S5 Table). Overall nattl2 had the highest expression in all morphs. The nattl1 had higher expression in embryos of PL-charr than in AC-and SB-charr, while nattl2 and nattl3 were more expressed in SB-embryos. Note however, the efficiency of the primers for the nattl genes ranged from 1.72 to 1.77, which suggests this data should be interpreted with caution. . 5 genes were differentially expressed between the two morphs (Alp, Krtap4-3, Lyz, Nattl, Parp6), while 4 further genes did not show significant expression differences between morphs (Cgat2, Cox6B1, Ndub6, Ubl5), see S4 Table. Error bars represent standard deviation calculated from two biological replicates.  Table. In order to evaluate the hypothesis that nattl genes have immunerelated functions we studied expression in adult tissues (in AC-charr). The nattl expression was highest in the gills, followed by expression in kidney, skin and spleen. Low expression levels were detected in liver, intestine and heart (S1 Figure and S5 Table). The three nattl paralogs followed different patterns, whilst each of them showed significant expression differences among tissues. Nattl1 was mainly expressed in spleen and kidney, while nattl2 showed a significantly higher expression in skin, liver and in gills. Similarly, the relative expression of nattl3 was highest in the gills and skin. This indicates that the three nattl paralogs are expressed in a tissue specific manner, and also differently during the development of the three charr morphs studied here.

Expression differences in the developing heads of benthic and limnetic charr morphs
To get a handle on the craniofacial divergence between sympatric Arctic charr morphs we used qPCR to study 8 paralog groups with expression difference in the RNA-seq data (all higher in SB). We focused on those with known craniofacial expression in zebrafish development 72 and compared two benthic (SB, LB) and two limnetic charr (AC, PL). We analyzed heads at three time-points (178, The eighth gene, Retinoic acid receptor gamma-A (Rarg) gave a small but significant response in the head, but the effects were reversed, i.e. the expression was higher in AC. The expression difference of the seven genes was, in almost all cases, consistent over the three timepoints studied (See S2 Figure). In summary the qPCR confirmed the differential expression of 12 of the 17 paralog groups studied (Table 3), some which had 5-10% FDR support. To us that suggests substantial expression differences between these two charr morphs, and that the data can lead to hypotheses about morph specific activity in particular structures, like the developing head.

Analyses of polymorphism in Arctic charr transcriptome
The RNA-seq data also revealed segregating variations with large frequency differences between charr morphs. To uncover candidate SNPs we mapped the reads to all of the S. salar EST-contigs. Filtering on coverage yielded 165,790 candidate SNPs (Table 4); of those 66.569 came from reads that mapped uniquely and 57.009 candidate SNPs from reads that mapped to more than one contig; with limited overlap between lists. Assuming that the expression of paralogous genes is stable, then differences among paralogs appear as SNPs at similar frequency in all samples. By requiring variant frequency differences (p < 0.05, uncorrected) between samples we reduced the list of candidates by two thirds, yielding over 20.000 candidate SNPs. Note, as cDNA from charr families was sequenced (not a population sample), estimates of SNP frequencies are imprecise. To err on the side of caution, we chose SNP candidates with 50% or higher frequency difference between morphs for further study. The candidate SNPs were also summarized by frequency of the derived allele, in reference to the S. salar sequence. This gave 672 and 872 SNPs at higher frequency, in AC-charr and SB-charr, respectively. The uniquely and multiply mapped reads, revealed approximately similar numbers of candidate SNPs. Gene ontology analysis showed that for derived SNPs in SB, there was an excess of variants in genes related to translation, both as a broad category and specific subgroups (S6 Table). There was also enrichment of SNPs in genes related to DNA-mediated transposition, DNA integration, DNA replication and oxidation-reduction process. No GO categories were enriched for high frequency derived SNPs in AC. Furthermore, functional effects of the candidate SNPs (UTR, synonymous and non-synonymous) were predicted. The distribution among those categories did not differ between variants detected by uniquely or repeatedly mapped reads, χ 2 [3] = 2.59, p = 0.46 (S7 Table).
A total of 60 candidate SNPs are nearly fixed in one morph, with frequency difference between morphs above 95% (after manual inspection of contigs and SNP position three candidates were removed since they represented the same SNP). Of these "fixed" SNPs 46 came from uniquely mapped reads and 14 from reads that mapped more than twice (Table 5 and Table 6). For the SNPs from uniquely mapped reads, 17 are fixed in AC-charr and 29 in SB-charr. The few genes with two or more polymorphic sites were; Keratin type II cytoskeletal 3 (Krt3), Cysteine sulfinic acid decarboxylase (Csad) and DNA-directed RNA polymerase I subunit RPA12 (Rpa12) with 5, 5 and 2 SNPs respectively. Krt3 and Csad had significant differentiation in both SB and AC. Similarly, 14 SNPs with large differentiation between morphs were predicted from reads that mapped on two or more contigs (Table 6). Of these, we found two variants in the mitochondrial 60S ribosomal protein L36 (RpL36) and variants in 4 other mitochondrial genes (28S ribosomal protein S18a mitochondrial (MRPS18A), Apoptosis-inducing factor 1 mitochondrial (AIFM1), Isocitrate dehydrogenase [NADP] mitochondrial (acIDH1) and Protein S100-A1 (S100a1)), all at higher frequency in AC-charr. PCR and Sanger sequencing of population samples confirmed SNPs in DNA2-like helicase (Dna2), a gene with nuclear and mitochondrial function, and two other genes Uroporphyrinogen decarboxylase (Urod), and Mid1-interacting protein 1-like (Mid1ip1) (S2 Table). The candidate variant Eukaryotic translation initiation factor 4 gamma 2 (Eif4g2) was not substantiated by the PCR/sequencing.

Polymorphism and expression of Arctic charr mtDNA
Considering the enrichment of differentially expressed genes related to mitochondrial energy metabolism (above), and high frequency Figure 6. Expression differences of craniofacial candidate genes in developing head of Arctic charr morphs. Relative expression ratios, calculated from the qPCR data, were subjected to an ANOVA to test the expression differences amongst four charr groups and three close time points (τs). The underlined gene names reflect significant difference between SB and AC-charr. A post hoc Tukey's test (HSD) was performed to determine the effects of morphs, time and morph-time interaction (M X T). White boxes represent low expression, while black boxes represent high expression. The shading represents significant different expression between the samples (α = 0.05, NS = not significant). The genes studied were, Claudin 4    candidate SNPs in several genes with mitochondrial function in AC-charr we decided to study the mitochondrial transcriptome further. The charr studied here reflect metabolic extremes, the aquaculture charr was bred for growth while the small benthic morph is thought to have experienced natural selection for slow metabolism and retarded growth 38,75 . Although mRNA preparation protocols were used for generating cDNA for the RNA-sequencing, a substantial number of reads came from non-polyadenylated sequences. By mapping the reads to mtDNA sequence of Arctic charr we could estimate expression and infer polymorphism both in genes and intergenic regions. There was a clear difference in sequencing coverage, with more than twice as many reads mapped from the AC-compared to SB-charr (mean fold difference 2.27, Wilcoxon test, p < 0.0004). Note, as only two types of fish are compared, the polarity of expression divergence is unknown.

(Cldn4), adseverin (Scin), Junction plakoglobin (Jup), Lipolysis stimulated lipoprotein receptor (Lsr), Major vault protein (Mvp), Transforming growth factor beta receptor II (Tgfbr2) Vitamin D receptor a (Vdra) and Retinoic acid receptor gamma-A (Rarg).
The mapped RNA-reads were used to identify polymorphism and divergence in the entire mitochondrial chromosome. The polymorphisms were found by mapping to mtDNA from a Canadian S. alpinus 48 , but ancestral vs. derived status inferred by comparison to S. salar mtDNA. This revealed 82 candidate sites, including 35 that represent divergence between Icelandic and Canadian charr. A total of 20 candidate SNPs had high (more than 50%) frequency difference between SB-and AC-charr ( Figure 7). There was no bias in the distribution of derived SNPs, 11 on the AC branch and 9 in SB. The divergence between Iceland and Canada is particularly little in the 12s and 16s ribosomal RNA genes. Curiously two SNPs in those genes differed strongly in frequency between morphs (Figure 7). To confirm and better estimate the frequency of variants in the ribosomal genes, we PCR amplified and sequenced two ~550 bp regions in the rRNA genes, comparing three morphs (PL, LB and SB) from Lake Thingvallavatn ( Figure 8A, C & E, S2 Table). The 12s polymorphism (m1829G>A) differed significantly between the morphs (χ 2 [2] = 8.6, p = 0.014), and was at highest frequency in the SB (0% in PL, 12.5% in LB and 75% in SB). Similarly m3411C>T in the 16s was enriched in SB (62.5%) but found at lower frequency in PL (0%) and LB (12.5%) (it differed significantly between morphs, χ 2 [2] = 9.3333, p = 0.009). The Sanger sequencing also revealed three other polymorphisms in the amplified region, not seen in the transcriptome. Among those m3211T>C in the 16s gene was at 75% frequency in LB, but not found in the other morphs (χ 2 [2] = 19.76, p < 0.0001).
In order to gauge the potential functionality of those variants we aligned the rRNA genes from nearly hundred fishes and several vertebrates. The position affected by m1829G>A and m3211T>C, in the 12s and 16s rRNAs, are not well conserved in fishes or vertebrates ( Figure 8B & D). However m3411C>T, in the 16s rRNA, alters a position that is nearly invariant in 100 fish genomes ( Figure 8F). The only exception is Pacific menhaden, which curiously also has T in this position. This region could not be aligned properly in other vertebrates. Thus m3411C>T alters a conserved position, but probably not very drastically as the introduced allele is tolerated in another fish species.  The file is tab-delimited and the columns are; "Unigene.Description": the annotation for that gene/paralog group. "NR.contigs": number of contigs with this annotation. "logCPM": count per million, log-scale. "logFC.morph": Mean fold change between the morphs, log-scale. "logFC.T163", "logFC.T200", "logFC.T433": Mean fold change for each timepoints compared to timepoint 141, log-scale. "FDR.morph": P-value for morph difference, multiple testing corrected. "FDR. time": P-value for time differences, multiple testing corrected. "Contigs": SalmonDB id for the contigs with the specific annotation 109 . "Gene Type": Designates the reference and candidate genes. "Gene": Name of the gene. "Morph": Which charr type the sample came from. "Relative age": Developmental timepoint, and also indicates the samples from adult fish. "Biological replicate": The two or more biological replicates used. "cDNA No": Marks the cDNA isolation used. "Ct value": Estimate of gene expression. "Sample": Indicates the material used, whole embryos or distinct tissues. "Batch": Demarcates distinct collections of cDNA, applies only to nattl 110 . "Gene Type": Designates the reference and candidate genes. "Gene": Name of the gene. "Morph": Which charr type the sample came from. "Relative age": Developmental timepoint. "Biological replicate": The two or more biological replicates used. "cDNA No": Marks the cDNA isolation used. "Ct value": Estimate of gene expression. "Tissue": Indicates the material used 111 .

Discussion
We are interested in the predictability of evolution at the molecular level, especially whether there exist principles that influence the rewiring of developmental and regulatory systems 4,76 . One way to study this is to identify genetic and developmental effects affecting key traits in species or populations which exhibit parallel evolution. The objective of this study were to get generate hypotheses about the genetic and molecular systems that associate with benthic morphology in charr by mainly focusing on the small benthic morph in Lake Thingvallavatn, Iceland. But as transcriptome were sequenced from embryos of SB-charr and aquaculture charr the data also reflect on the genetics of charr domestication.
Developmental transcriptome of Arctic charr morphs As no reference genome is available for Arctic charr, we mapped reads to S. salar EST-contigs 57 in order to estimate expression and identify candidate genetic polymorphisms. As many of the contigs are short or have overlapping annotations, we collapsed genes into paralogous genes when appropriate for the expression analysis. The main advantage of this approach was the reduction of the number of statistical tests (and hence an increase in statistical power). The downside is that paralog-specific expression patterns are masked, as our qPCR results of the natterin like gene family show ( Figure 5 and S1 Figure). Recent rainbow trout data shows about 1/4 of paralogs from the latest whole genome duplication event retain the very similar expression patterns 22 indicating that distinct expression patterns of two paralogs is quite common 77 . In their analysis of the Arctic charr gill transcriptome, Norman et al. (2014) 13,14 also used Illumina sequencing technology to evaluate expression. Their reads were longer (2x100 bp) than in this study (36 bp) enabling them to assemble contigs. They did not consider the paralogs in their approach and merged contigs based on sequence identity. Thus the complexity of Arctic charr transcriptome still remains unsolved. Our data reflects differential deployment of several gene classes during Arctic charr development. Studies in salmonids and other fish have demonstrated large changes in expression during early development, including coordinated changes in many cellular and developmental systems 9,78-81 . Several blood coagulation factors genes showed significant changes during charr development, and were also more highly expressed in the SB-charr. This might reflect differences in the rate of development of blood composition, or tissue composition, in the two morphs. While our main interest is on the derived and repeatedly evolved small benthic charr, the data can also reflect differences due to domestication. In this study we chose to compare SB to AC-charr for several reasons, i) AC-charr has limnetic like head morphology, ii) was available for harvesting of running fish, and iii) because we wanted a strong contrast in this first survey of charr developmental diversity. The AC-charr proved a useful, as the data presented here has already revealed differential expression of several developmental genes and regulators with differential expression between benthic and limnetic charr 51,52 . Previously we found tight correlation of RNA-seq expression and qPCR estimates -using data from this very transcriptome 51 . Furthermore, we have actually used the same morphs (AC and SB) and samples in a comparison of the developmental miRNA transcriptomewhich reveal that expression of several miRNAs correlates with morph differences 56 .
Higher expression of lysozyme II C and natterin-like in SB-charr Natural selection can shape variation in immunological genes. We decided to study further Lyz2 and the putative immunological genes nattl that had higher expression in SB. The substrate of lysozyme 82 is the bacterial cell wall peptidoglycan and it acts directly on Gram-positive bacteria 83 . Lysozyme also promotes the degradation of the outer membrane and therefore indirectly acts also on Gram-negative bacteria 84 . Another gene that caught our attention was natterin-like. Natterins were first discovered from the venom gland of the tropical toxic fish species Thalassophryne nattereri 70,71 , and are found by sequence similarity in e.g. zebrafish, Atlantic salmon and here in Arctic charr. The Natterin proteins contain a mannose-binding lectin-like domain (Jacalin-domain). Mannosebinding lectins are pathogen recognition proteins (antibodies) and therefore are important for the acute phase response of fish 85,86 , thus we hypothesized that nattl genes in charr may have immune related functions. The data are consistent with this as the highest expression was found in skin and kidney. This putative immune functions needs to be verified. It is possible that higher expression of those two genes in SB-charr reflect preparation of juveniles for bottom dwelling habitats, which may be rich in bacteria and challenging for immune systems. One can ask whether immunological genes are expected to show similar or less parallelism than others genes shaped by natural selection?
The current data does not reflect on this question, but our population genetic work shows genetic variation in immunological genes (MHCIIα and cath2) does not correlate with the SB-charr ecotype in Iceland 45 .
In this study we collapsed contigs into paralog groups for the transcriptome analyses. The disadvantage of this approach is that differential expression of a paralog, can be masked by related genes that do not differ between groups. We looked at this by studying the expression of three paralogs of the natterin like genes in different morphs during Arctic charr development, and among tissues of adult AC-charr. The data suggest that the three nattl genes are expressed differentially between the morphs, thus it is not divergence in the expression of one paralog that explains the general nattl expression disparity in the transcriptome. Certainly, other scenarios could apply to other genes in the transcriptome.

Expression divergence in craniofacial genes in benthic morphs
A study of the skulls of charr post-hatching embryos and juveniles from Lake Thingvallavatn, showed that some elements of the developing head ossified earlier in SB-charr than in PL-charr 87 . Morphometric analyses of developing heads (same stages as studied here) demonstrate differences in craniofacial elements between AC-and SB-charr, along a limnetic vs. benthic axis 74 . Based on those developmental phenotypes we investigated further genes with roles in craniofacial development that were differentially expressed in the transcriptome. Guided by this transcriptome we had already found two extra-cellular matrix (ECM) remodeling genes, Mmp2 and Sparc and a conserved co-expression module of genes with known roles in craniofacial morphogenesis, to have higher expression in developing heads of benthic Arctic charr morphs than in limnetic morphs 51,52 . Bioinformatic and qPCR analyses suggest the co-expression module may potentially be affected by quantity of the transcription factor ETS2. These studies and the current data confirm the utility of the contrasting developmental transcriptomes for identifying candidate genes with differential expression during head development, as 7 out of 8 candidates were confirmed by qPCR. These genes had consistently higher expression in the developing head of two benthic morphs (SB and LB), and lower in more limnetic fish (AC and PL). The most noteworthy aspect is the fact that three of the morphs (SB, LB and PL) are closely related and live in sympatry in Lake Thingvallavatn 44 .
We focused on a few targets of Tgf-β and Ahr signaling pathways because of their role in craniofacial morphogenesis and transcriptional connection 88-90 . Adseverin (Scin) was one of the top differentially expressed genes (Table 1) and has roles in rearrangements of the actin cytoskeleton, chondrocyte differentiation and skeletal formation 91,92 . Also, in the transcriptome Lsr, Cldn4 and Tgfbr2 had higher expression in SB-charr, and we show that higher expression of those genes associated with the benthic morphotype. Lsr is a molecular component of tri-cellular tight junctions 93 and has been shown to be suppressed upon Tgf-β1 stimulation 94 in a human cell line. Similarly, Cldn4, a tight junction protein with unknown role during embryonic morphogenesis, is a target of the Tgf-β and Ahr signaling pathways 95,96 . Finally, the expression of Tgfbr2, encoding a receptor of Tgf-β was slightly but significantly higher in the head of benthic morphs. Previous studies suggest a crucial role of Tgfbr2 in craniofacial morphogenesis 97 .
We also confirmed differential expression of other genes, including two with higher expression in SB-charr. Mvp is the predominant component of cytoplasmic ribonucleoprotein structures called vaults 98 , which is highly conserved across eukaryotes. The vaults have been something of an enigma, but are implicated in several processes from signal transmission and immune response 99 . The Jup gene also showed higher expression in SB-charr. Finally, higher expression of Vdra, encoding the vitamin D receptor A, was found in the heads of benthic charr. The receptor regulates mineral homeostasis, osteoblast differentiation and bone metabolism 100 . A related study from our group, also building on this dataset, mapped in more detail the differential expression of these and other coexpressed genes in limnetic and benthic charr 53 .
To summarize, the results show that RNA-sequencing of Aquaculture charr with limnetic craniofacial morphology and small benthic charr can be used to reveal differential expression of genes that associate with limnetic and benthic divergence in craniofacial elements in sympatric charr morphs. It would be interesting if expression of these genes associates with benthic morphology in independently evolved charr populations, as was seen for certain mTOR-pathway genes in muscle of adult SB-charr 47 , or even in other species with similar trophic diversity.
Genetics differences between the AC and SB-morphspossibly in mtDNA function Previous studies on microsatellite markers documented the population history of charr in Iceland and in particular the parallel evolution of SB-charr 44 . Our data confirm genetic differences between SB and AC-charr. By comparing AC and SB-charr, that represents a small benthic resource morph that has evolved repeatedly in Icelandic stream and pond habitats 44 , we hoped to implicate genes and pathways involved in adaptation to these special habitats. But the AC-charr is also interesting, as domestication over several decades has led to rapid growth and increased size 50 . Morphometrics have not been used to compare the body or craniofacial shape of AC to other charr morphs, but domestication of O. mykiss has affected body shape and fin structure in partiuclar 101 . The allele frequency differences and expression divergence observed can reflect neutral population genetic processes and/or selection during domestication or adaptation of SB-charr. By studying expression and allele frequencies in limnetic and benthic morphs from more locations, it may be possible to disentangle these questions. We restricted ourselves to verification of several SNPs, and focused mostly on variants in mtDNA because to us the data suggest interesting divergence in systems related to energy metabolism. First, there is 2X higher expression of respiratory electron transport chain components in AC compared to SB-charr and 100% more mitochondrial derived reads are found in the AC-charr samples. Note that the direction of divergence is unknown, i.e. whether expression was up in AC or down in SB. Second, many derived candidate-SNPs in genes related to mitochondrial function were at high frequency on the AC branch. For instance in S100A1, which has been implicated in mitochondrial regulation in cardiac tissue in humans 102 , but its expression is probably not exclusive to this tissue. Third, while the mitochondrial ribosomal genes generally evolve slowly, we do see derived variants at high frequency in the SB and large benthic charr in Lake Thingvallavatn. Specifically, m3411C>T in SB affects a position that is highly conserved among fish, and could affect function of the 16s rRNA. Earlier studies of mitochondrial markers in S. alpinus did not find large signals of divergence within Iceland 40,42,45 , probably because they studied other genes.
The mitochondrion is more than a powerhouse, it integrates metabolism, cell cycle and apoptosis 103 . The number of mitochondria and its functions are known to correlate with environmental attributes. For instance in Antarctic fishes under extreme cold, higher numbers of mitochondria are found in muscle and heart cells 104 . Our data suggest an expression difference between morphs that could reflect differences in total number of mitochondrion, the number of mtDNA copies per mitochondrion or cell, or difference in RNA expression from the mtDNA, possibly due to evolution of mtDNA related to diet and/or temperature 105 . In sum, the results suggest divergence (adaptive or neutral) in mitochondrial function, due to the domestication of aquaculture charr and/or adaptation of the small benthic charr to its habitat in Lake Thingvallavatn. But further work is needed to map out the expression differences of mitochondrial related genes in more SB and anadromous charr morphs (representing the ancestral state). The mtDNA signals could also be investigated in populations along ecological clines (e.g. temperature) or with respect to life history 106 .

Conclusions
The data presented here suggest genetic and expression changes in multiple systems associate with divergence among the highly polymorphic and rapidly evolving Arctic charr in Iceland. The data reveal differential expression of two immunological genes between morphs and of several craniofacial developmental genes, that may help sculpture benthic vs. limnetic heads. The genetic data suggest among other things differentiation in the charr mtDNA between the SB and AC-charr morphs. It must be acknowledged that it is not trivial to identify genes affecting variation in ecologically important phenotypes, like shape 107,108 . Our broad interest is in how natural selection tweaks genetic regulatory systems, for instance via genetic changes in regulatory sequences or post transcriptional modifiers relating to adaptations. Genetic changes affecting gene expression can be raw material for adaptation, but could also rise in frequency due to reverberations in regulatory cascades 76 . Following this work we plan to study the degree of developmental and population genetics parallelism of the small benthic charr, typically found in cold springs and small pond habitats in Iceland with lava substratum 29,44 . The availability of charr populations at different stages of divergence sets the stage for future genomic studies of the roles of genes, environment and plasticity for shaping this polymorphic species.
• Analyses of RNA sequencing data: JG, AP.

Competing interests
No competing interests were disclosed. I confirm that the funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant information
Supporting Information S1 Figure.
S2 Figure. Relative expression of selected craniofacial candidate genes. Relative expression of 12 candidate genes with characterized craniofacial expression during zebrafish development (ZFIN website) in the head of SB, LB, PL and AC at three time points in development.
In the transcriptome data all of the genes had shown higher expression in SB at 200 τs. The expression is normalized to the geometric means of two craniofacial reference genes (ACTB and IF5A1). Expression is relative to a replicate of AC morph at 200 (τs), set to one. Error bars represent standard deviation calculated from two biological replicates and each biological replicate contains homogenate of six heads. From RNA-reads that mapped to one (Uni) or more (Rep) S. salar ESTs.

Supplemental
The candidate SNPs frequencies differ more than 50% between SB and AC-charr, summarized by which morph with higher frequency of the derived allele. The study by Gudbrandsson et al. reports a thorough analysis of differences in the transcriptome between different 'morphs' or 'populations' of artic charr. More specifically they have studied the transcriptome of eggs and larvae from a natural population of small benthic charr (SB) and Icelandic aquaculture charr (AC), which is fast growing and have a 'limnic-like' morphology. They find a list of potential candidate genes involved in the ecological differentiation of artic charr (and during the embryonic development). In addition they studied the transcriptome from different tissues of adult AC-charr. From the transcriptome of these populations two populations they developed 12 SNP-markers applied to other sympathric (Lake Thingvallavatn) wild morphs to study if these genes differed between other morphs. Finally they also study mtDNA expression between morphs to find they mainly differ between the benthic morphs and a limnic.
The search for genes involved in the ecological divergence of species is an important topic that has exploded the last decade with the new generation of sequencing. I find this study to be an important contribution because of the study system with artic charr is an example of relative recent and rapid divergence into many different morphs/ecological, and the extensive and thorough investigation of the differences in the transcriptome between morphs.
However, I think the authors try to stretch their conclusions a bit too far. The study is great as a base for further research in the topic, which I guess is in the pipeline. The use of cultivated charr make sense for comparing the most extreme morphs. But to me it does not make sense for making conclusions about genes involved in the ecological niche differentiation in natural populations, which is the motivation of the study in the introduction and brought up in the discussion). The cultivated population has been selected for fast body growth and they are not from Lake Thingvallavatn and little can therefore be said about the genetics of the ecological differentiation of sympatric species. Not very surprising genes related to metabolism seemed upregulated in AC and immunogens upregulated in SB. What does that actually tells us about the genetics of ecological differentiation of natural populations??
Although the aims on p. 4 feels valid, they are not contingent with the previous text in the introduction. Thus, I suggest that the much of the earlier part of the introduction is rewritten to actually address the differences in gene expression between a cultivated morph and its extreme opposite small benthic artic charr.
As far as I understand it is only egg that are kept in the same environment, but the parents have been raised in different environments and transgenerational plasticity cannot be ruled out. This is not a major criticism (the ideal case would be to have had them in lines in a common environment of course) but F1000Research raised in different environments and transgenerational plasticity cannot be ruled out. This is not a major criticism (the ideal case would be to have had them in lines in a common environment of course) but needs to be addressed in the text.
The 'Nattl' paralogs provide an interesting case where the expression of different paralogs has been studied. But again, are the result difficult to interpret from an ecological niche differentiation perspective. Often is the natural small limnic morph (PL) in between SB and AC (Fig. 5A), but for Nattl1 and Nattl2 AC and SB seem most similar? The connection to the original question is weak and the authors do not conclude more than "…it is not divergence in the expression of one paralog that explains the general nattl expression disparity in the transcriptome." Fair enough, but that is more about the genetic architecture than the genetics of ecological divergence.
In the validation of the transcriptome differences with qPCR of 9 genes/paralogs only 5 was still significant. What conclusions should one make out of that, that around half of the 296 paralogs differing between SB and AC are false detection (despite FDR < 5%). I support the use of qPCR but please comment on the implication of this.
I think Fig. 6 should be converted into a bar-graph plot instead (this feels more like a table).
To conclude, I think this study is a great contribution for suggesting putative differences in the transcriptome of a fish species. However, the importance for understanding ecological differentiation of sympatric species it is, however, so far limited as it that would require more using natural morphs, replicated populations, back-crosses, investigation of plasticity and how reproductive isolation is maintained etc., which is likely to come. But until that, I suggest the this text should mainly focus on the differences between SB and AC arctic charr, and not try to squeeze in everything in one paper. Note that on p. 8 the qPCR of 8 paralogs in embryonic heads is metioned but the results do not come until "Expression differences in the developing heads of benthic and limnetic charr morphs".

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
No competing interests were disclosed.  Introduction: " " should be " In this study, we compare SB-charr from In this study, we compared " (again, it is correct here to use past tense -the authors should check the rest of the SB-charr from manuscript for similar tense issues). . have tested for differential gene expression at multiple developmental time-points et al among a number of Artic charr morpho-types from Lake Thingvallavatn (3 wild morphs, 1 studied with RNA-seq and qPCR, the others with qPCR only) and Holar aquaculture (1 domesticated morph, RNA-seq and qPCR). They have also studied multiple tissues/body regions for a subset of the differentially expressed genes found with RNA-seq. The goal of the paper was to find candidate genes that may underlie variation in morphology, with a focus on craniofacial morphology related to benthic vs. limnetic feeding. In general, I think this goal was met and this paper contributes to our understanding of the mechanisms contributing to morphological evolution in a non-genetic model organism. The authors provide an extensive, multi-time point comparison of two morphologically divergent groups of charr reared in a common environment (reducing the influence of phenotypic plasticity) and have collected a 1.
provide an extensive, multi-time point comparison of two morphologically divergent groups of charr reared in a common environment (reducing the influence of phenotypic plasticity) and have collected a tremendous amount of data. This information will help them to hone in on the genetic loci contributing to phenotypic evolution in this very interesting system, and on the effects of domestication. However, there are a number of major issues that do need to be more clearly addressed in the manuscript prior to final publication. I have outlined these comments below.

Major Comments
Introduction: Requires some reorganization, clarification of what phenotypes have evolved in parallel among morphs, and how the authors separate the effects of domestication (SB vs. AC) from benthic/limnetic evolution (SB/LB vs. PL/AC). a) At present, the introduction focuses upon the utility of instances of parallel evolution to help us determine how repeatable evolutionary change may be. This is definitely true, and the repeated evolution of the dwarf, benthic morph (SB; the focus of the introduction/abstract/discussion) in many lakes strongly argues that this phenotype has evolved via natural selection. However, it is not clear to me if true 'parallelism' seen among the SB (small benthic) and LB (large benthivorous) vs. AC (Holar aquaculture) and PL (small planktivorous) morphs because not enough information is provided for me to assess this. To support the argument for parallelism the specific traits that have evolved in parallel among morphs must be displayed and the evolutionary history of these morphs should be clarified (e.g. in paragraph 6 and Figure 1). As well, any related non-parallelism in traits should also be discussed (i.e. how are the domesticated AC and wild PL different?). At present Figure 1 only shows the AC and SB morphs, and does not point out the specific traits they are interested in. This is critical background information for readers who are not familiar with this system.
b) The comparison of AC (domestic, limnetic-like head) vs. LB (wild, benthic like head) looks at two confounded variables: domestication and the benthic/limnetic morphology. This should be clearly stated in the introduction, and the use of the additional morphs (PL, LB) in detangling domestication vs. benthic/limnetic evolution should be noted.
c) The use of the AC morph is still a bit unclear to me. The argument for point 'ii) of the availability of abundant AC material' could be expanded by providing more information on the 'limnetic' like features of this morph and why it is an appropriate comparison to a benthic morph, the genetic divergence from the lake Thingvallavatn fish, and also the selection regime it has experienced (selection for limnetic features? What other traits vary with domestication?). d) Paragraph 2 -Much of this paragraph, including discussing the ability to measure gene expression and relate to phenotype in fishes, is unnecessary as fish are no different from other vertebrates in this respect. Instead, the final sentence "One approach to identify pathways related to function or morphological differences is to study gene expression during development" should become the 'topic sentence' and expanded upon to explain why gene expression studies are especially relevant ways to link genotype to phenotype in evo-devo studies. e) Better highlight the strengths -The authors have done a wonderful job of assessing multiple developmental time points and rearing fish in a common garden environment. However, they do not highlight these strengths. Some small notes on the importance of controlling for phenotypic plasticity in these traits (which are known to be quite plastic) to better study genetic differentiation would be a nice addition.
would be a nice addition.
Methods: a) Page 4 paragraph 1 -Clarify the number of fish used to make the crosses (this will help us determine the likelihood of selecting a full or half-sib for sequencing/qPCR).
b) I should note that I am not an expert in the analysis of RNA-seq data, but luckily the first reviewer has done an excellent job of commenting upon these aspects of the project. I fully agree with their comments and suggestions. I would also like to see more information on the methods used to pool samples and how RNA-seq data was normalized among samples, developmental times and morphs. I will also note that the authors often use for comparisions, not , which is S.salar O.mykiss a closer relative to . The reasons for this approach should be discussed. S.alpinus c) I am also not trained as a population geneticist. However, from my experience studying paralogous genes in salmonids, and with respect to the author's own findings for the Nattl paralogs (Fig 4), I do not think it is prudent to "assume that the expression of paralogous genes is stable… " in the methods (page 12). In fact, Berthelot . (2014)  d) The authors should use their genetic information to test if the fish chosen are siblings with each other (full or half-sibs). This may have important implications for the population genetic analyses. e) Page 5 -It is not appropriate to change the meaning of the word 'gene'. I think it is much clearer to use the term 'paralog group' or 'gene family' when referring to the fact that the authors do not study single genes, but instead groups of paralogs. f) Selection of genes for qPCR -the methods by which genes for the qPCR studies (Fig 3) were selected should be clearly noted. From my reading, it seems that most of these genes do not significantly vary among SB and AC at the 1% FDR level (Tables 1 and 2; only Natterin?). Thus, I am assuming these genes are only significant at the 5% FDR level (S1 file) -why focus upon these and not those significant at 1%? As well, it would be good to include information on why different genes were selected for Figure 3 (qPCR validation of whole fish) and Figure 4 (candidate genes-qPCR validation in just the head). Finally, the abbreviations used for qPCR validation should also be listed in Table 1 for easy comparisons among figures/tables.

Results & Figures:
a) Include an experimental design figure -At present, it is difficult to keep track of all of the morphotypes, tissues, and developmental time points used without referring to the methods. Thus, an experimental design figure summarizing the samples used (morphotype, population, sample size, developmental time point), how they were pooled and which techniques were used to measure gene expression on each sample (RNA-seq and/or qPCR) is needed. b) Include the LB and PL morphs in Figure 1 and clarify traits of interest -The legend states that "differences in size, coloration and head morphology are apparent", but it would be better to specifically point out the differences they are referring to. F1000 is for a general audience, and this would help non-ichthylogists better understand what ecologically-important traits the authors are interested in (e.g. those related to benthic/limnetic feeding). In addition, the two other morphs used in the qPCR studies should also be displayed (large benthivorous and small planktivorous) to in the qPCR studies should also be displayed (large benthivorous and small planktivorous) to facilitate phenotypic comparisons and assess parallelism in benthic/limnetic feeding and/or the effects of domestication on AC. c) Figure 5-this is actually a table not a figure (?) and is a bit confusing. I think it is much easier to interpret Figure S2 (displaying the data as in Fig 3 and 4), and that Fig 5 and S2 should be switched. It would be great to show significant differences in mRNA content in this, and all other figures, by including symbols. Also, full gene names should be listed in all figure legends.

Discussion:
The discussion focuses on the SB morph (page 17 -"The objective of this study were to get a handle on genetic and molecular systems that associate with benthic morphology in charr by mainly focusing on the small benthic morph in Lake Thingvallavatn, Iceland"), while the introduction discusses parallel evolution (indicating that the comparisons should be among many morphs). These are two different topics i) mRNA content differences among benthic vs. limnetic morphs changing in parallel or ii) linking mRNA content to phenotype in SB (benthic, wild) vs. AC (limnetic head, domesticated) morphs. In particular, the role of domestication vs. wild fish divergence needs to be addressed. At present these two topics/questions are mixed in the introduction/discussion and should be addressed separately. a) Paragraph on Immune Defences -Is immunity also expected to evolve in parallel in all benthic morphs? Is this predicted to be unique to SB vs. AC? Whatever the case, the parallelism (or not) in these genes should also be discussed, and whether this relates more to domestication in AC or differences between limnetic vs. benthic fish. Much of the functional discussion can also be cut. b) Page 18 -The information about genes found to be differentially expressed among morphs in your prior work should also be in the introduction, as it is background work that explains why you took this transcriptomic approach. This can also be used to explain why you focused in on particular qPCR genes. c) A discussion of domestication related differences vs. benthic/limnetic differences should be included. I think the data from head gene expression is very interesting (Figs 5, S2) and really speaks to this question. d) In general, the role of stochastic evolutionary processes, and not just selection (artificial and natural) should be noted. For example, if the AC charr were simply taken from a stock with a different mtDNA haplotype then these differences in the mtDNA genome might not be adaptive, just random. If the AC fish has much higher mtDNA expression might this be simply a domestication issue and not indicative of selection in SB as stated? Finally, you find that not all mitochondrial transcripts (which are transcribed as a polycistronic transcript) are found at similar levels (Table 1) -what does this tell you about differential degradation/post-transcriptional processes? e) There is no discussion about the "Analyses of polymorphism in Arctic charr transcriptome" (Table 3, 4, 5), except for the mtDNA.

Minor Comments
Introduction: F1000Research Introduction: a) Paragraph 3 -"Furthermore, recent estimates from the rainbow trout….by utilizing multiple data sources the genome assembly problem of this family can be solved". I am not sure how this statement is relevant to this particular study. This and the following statement seem more appropriate for the methods/discussion to me.
b) The morphs being discussed should be clarified throughout the paper. For example, the authors often state "among morphs/among charr populations" but it is not clear which of the many morphs they are referring to (e.g. Paragraph 5, first sentence on allozymes and mtDNA and later sentence on MCHIIa -do you mean all 4 morphs of specific 2-way comparisons? Are some morphs more differentiated than others?) Methods: a) The authors should note why they did not use the PI (large piscivorous) morph in any qPCR studies (in the methods or discussion) as this would be a nice morph to use in their tests for parallelism. b) Page 5 (last paragraph) -the methods used to remove particular variants needs to be clarified. In particular, why the assumptions used to remove variants are valid by referencing past studies. c) Figure 7 -It is not clear to me which variant is present in which morph. Adding the nucleotide to the x-axis (i.e. frequency of m1829G for B) would make this figure easier to quickly interpret. The "A.charr_WT" and "A.charr_M" should also be defined in the legend and it would be more appropriate to use scientific names for all species.
Discussion: a) Discussion of reference 32 -The discussion of reference 32 is not put into the proper context. Figure 6 of this paper (Berthelot . 2014) shows that there are many genes that have no correlation among et al expression patterns and/or differences in expression levels (1573, 1248, and 1895=4716 paralog pairs), and that together these represent more than the 1,407 correlated/similar expression level paralogs. This section of the discussion needs to be modified.
b) The Norman . (2014) paper should be mentioned earlier -if this is available why was it not used for et al their analyses? As well, the last sentence in this paragraph can be cut as it is evident. c) Page 18 -"Our new data also demonstrate differences in craniofacial elements between AC-and SB-charr, along a limnetic vs. benthic axis ". Are you referring to ref 79 or data from this study? If you are referring to 79, clarify and note what you found. This occurs a few times in the discussion General grammatical errors 79 F1000Research

General grammatical errors
There are a number of grammatical errors throughout this paper (e.g. "31 genes were higher expressed in SB and 40 genes higher in AC-charr"; "that may help sculpture benthic vs. limnetic heads" pg 19).
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
No competing interests were disclosed.

Competing Interests:
Author Response 04 Apr 2016 , University of Iceland, Iceland Arnar Palsson

Major Comments
Introduction:

Requires some reorganization, clarification of what phenotypes have evolved in parallel among morphs, and how the authors separate the effects of domestication (SB vs. AC) from benthic/limnetic evolution (SB/LB vs. PL/AC).
a) At present, the introduction focuses upon the utility of instances of parallel evolution to help us determine how repeatable evolutionary change may be. This is definitely true, and the repeated evolution of the dwarf, benthic morph (SB; the focus of the introduction/abstract/discussion) in many lakes strongly argues that this phenotype has evolved via natural selection. However, it is not clear to me if true 'parallelism' seen among the SB (small benthic) and LB (large benthivorous) vs. AC (Holar aquaculture) and PL (small planktivorous) morphs because not enough information is provided for me to assess this. To support the argument for parallelism the specific traits that have evolved in parallel among morphs must be displayed and the evolutionary history of these morphs should be clarified (e.g. in paragraph 6 and Figure 1). As well, any related non-parallelism in traits should also be discussed (i.e. how are the domesticated AC and wild PL different?). At present Figure 1 only shows the AC and SB morphs, and does not point out the specific traits they are interested in. This is critical background information for readers who are not familiar with this system.
These are excellent suggestions. At the end of the intro we stress the difference between Reply: the aims of our research program (study the genetics of parallel evolution) and the aims of this study (get a handle on differences between sympatric morphs, with the AC as possible outgroup). The morphs studied here do not represent parallel evolution of benthic phenotypes (SB and LB are both from the same lake and appear to be closely related -Kapralova et al 2011). Analyses of that question requires further studies. This data can implicate genes that separate PL/AC and SB/LB and may be studied in such follow up analyses of more populations. We have updated figure 1 as advised -including the 4 morphs studied, expanded on the legend and also provide an overview of research approach (part B).
b) The comparison of AC (domestic, limnetic-like head) vs. LB (wild, benthic like head) looks at two confounded variables: domestication and the benthic/limnetic morphology. This should be clearly stated in the introduction, and the use of the additional morphs (PL, LB) in detangling domestication vs. benthic/limnetic evolution should be noted.

F1000Research domestication vs. benthic/limnetic evolution should be noted.
c) The use of the AC morph is still a bit unclear to me. The argument for point 'ii) of the availability of abundant AC material' could be expanded by providing more information on the 'limnetic' like features of this morph and why it is an appropriate comparison to a benthic morph, the genetic divergence from the lake Thingvallavatn fish, and also the selection regime it has experienced (selection for limnetic features? What other traits vary with domestication?).
(b and c): The reviewer is correct, AC and SB are separated by multiple traits, and the data Reply probably reveal signals associating with most of them. Unfortunately the AC charr is not well characterized phenotypically, thus we can not address the question of other traits. We focus mainly on the head and jaw morphology, as these attributes distinguish benthic and limnetic morphs. The revised intro elaborates on the choice of AC, and how the follow up work on the morphs from Lake Thingvallavatn can help us sort this out. This point is also picked up in the discussion. d) Paragraph 2 -Much of this paragraph, including discussing the ability to measure gene expression and relate to phenotype in fishes, is unnecessary as fish are no different from other vertebrates in this respect. Instead, the final sentence "One approach to identify pathways related to function or morphological differences is to study gene expression during development" should become the 'topic sentence' and expanded upon to explain why gene expression studies are especially relevant ways to link genotype to phenotype in evo-devo studies.
: We restructured and shortened this paragraph around this topic sentence -and gave more Reply room for the previous RNAseq study on Arctic charr. e) Better highlight the strengths -The authors have done a wonderful job of assessing multiple developmental time points and rearing fish in a common garden environment. However, they do not highlight these strengths. Some small notes on the importance of controlling for phenotypic plasticity in these traits (which are known to be quite plastic) to better study genetic differentiation would be a nice addition.
: Great advice, we tried to integrate this into the last paragraph of the intro. Reply

Methods:
a) Page 4 paragraph 1 -Clarify the number of fish used to make the crosses (this will help us determine the likelihood of selecting a full or half-sib for sequencing/qPCR).
: We did bulk crosses, joining eggs from 5-10 females in a can and sperm from 3-5 males Reply (SB, PL, LB) and single parent cross for AC. Each sample included RNA pooled from 3 embryos, so there is a chance that full sibs were sequenced, but unlikely. The embryos/samples for qPCR are from similar pools. Now described better in methods.
b) I should note that I am not an expert in the analysis of RNA-seq data, but luckily the first reviewer has done an excellent job of commenting upon these aspects of the project. I fully agree with their comments and suggestions. I would also like to see more information on the methods used to pool samples and how RNA-seq data was normalized among samples, developmental times and morphs. I will also note that the authors often use S.salar for comparisions, not O.mykiss, which is a closer relative to S.alpinus. The reasons for this approach should be discussed.

F1000Research
: The RNA was isolated from individual embryos, quantified and then united (in equal Reply concentrations) prior to cDNA synthesis. The read counts per gene are normalized per million reads in sample. Not normalized with other variables. c) I am also not trained as a population geneticist. However, from my experience studying paralogous genes in salmonids, and with respect to the author's own findings for the Nattl paralogs (Fig 4), I do not think it is prudent to "assume that the expression of paralogous genes is stable… " in the methods (page 12). In fact, Berthelot et al. (2014) find the opposite (see my comments for the discussion).
: Excellent suggestion. We corrected our misunderstanding, added this fact into the intro and Reply discussion, and reinterpreted our data in this light.

d) The authors should use their genetic information to test if the fish chosen are siblings with each other (full or half-sibs). This may have important implications for the population genetic analyses.
:The fish chosen for pop-gen work are random sample from spawning grounds -assumed to Reply be not sibling groups. Our earlier study (Kapralova 2011) showed no family structure in charr collected this way from the lake.

e) Page 5 -It is not appropriate to change the meaning of the word 'gene'. I think it is much clearer to use the term 'paralog group' or 'gene family' when referring to the fact that the authors do not study single genes, but instead groups of paralogs.
: Excellent suggestion. We amended this., and use paralog group throughout. Reply f) Selection of genes for qPCR -the methods by which genes for the qPCR studies (Fig 3) were selected should be clearly noted. From my reading, it seems that most of these genes do not significantly vary among SB and AC at the 1% FDR level (Tables 1 and 2; only Natterin?). Thus, I am assuming these genes are only significant at the 5% FDR level (S1 file) -why focus upon these and not those significant at 1%? As well, it would be good to include information on why different genes were selected for Figure 3 (qPCR validation of whole fish) and Figure 4 (candidate genes-qPCR validation in just the head). Finally, the abbreviations used for qPCR validation should also be listed in Table 1 for easy comparisons among figures/tables.
: Very important point. We deliberately studied some genes with less statistical support (FDR Reply between 5% and 10%), to gauge the differences in the genes with less support and in particular to have a bigger pool of candidates that may relate to the specific developmental process (like head and jaw formation). Of course we can not assert that all the genes with strongest DE signal in the transcriptome are true positives, but the data can be used for hypothesis generation. We also amended table 1 and the figure legends accordingly. F1000Research measure gene expression on each sample (RNA-seq and/or qPCR) is needed. Figure 1 and clarify traits of interest -The legend states that "differences in size, coloration and head morphology are apparent", but it would be better to specifically point out the differences they are referring to. F1000 is for a general audience, and this would help non-ichthylogists better understand what ecologically-important traits the authors are interested in (e.g. those related to benthic/limnetic feeding). In addition, the two other morphs used in the qPCR studies should also be displayed (large benthivorous and small planktivorous) to facilitate phenotypic comparisons and assess parallelism in benthic/limnetic feeding and/or the effects of domestication on AC.

b) Include the LB and PL morphs in
: (a and b) Excellent suggestions. Now picture 1 has all 4 morphs, and a schematic Reply describing the work flow and samples. c) Figure 5-this is actually a table not a figure (?) and is a bit confusing. I think it is much easier to interpret Figure S2 (displaying the data as in Fig 3 and 4), and that Fig 5 and S2 should be switched. It would be great to show significant differences in mRNA content in this, and all other figures, by including symbols. Also, full gene names should be listed in all figure legends.
: We acknowledge that this graph is not the simplest, but would like to keep it over Figure S2.

Reply
Our reasoning is that this graph illustrates the sharp differences between the limnetic (AC-PL) and benthic (SB-LB), which are the main result in this section. But we will of course switch them, or possibly join both in a single figure ?? if the reviewer insists or the editors recommend it.

Discussion:
The discussion focuses on the SB morph (page 17 -"The objective of this study were to get a handle on genetic and molecular systems that associate with benthic morphology in charr by mainly focusing on the small benthic morph in Lake Thingvallavatn, Iceland"), while the introduction discusses parallel evolution (indicating that the comparisons should be among many morphs). These are two different topics i) mRNA content differences among benthic vs. limnetic morphs changing in parallel or ii) linking mRNA content to phenotype in SB (benthic, wild) vs. AC (limnetic head, domesticated) morphs. In particular, the role of domestication vs. wild fish divergence needs to be addressed. At present these two topics/questions are mixed in the introduction/discussion and should be addressed separately.
: We tried to separate these two aims more clearly in the revised discussion. The strategy Reply was to use the AC vs SB contrast for hypothesis generation, as the first aim is central to our program. We have now added sentences on the domestication in two parts of the discussion. a) Paragraph on Immune Defenses -Is immunity also expected to evolve in parallel in all benthic morphs? Is this predicted to be unique to SB vs. AC? Whatever the case, the parallelism (or not) in these genes should also be discussed, and whether this relates more to domestication in AC or differences between limnetic vs. benthic fish. Much of the functional discussion can also be cut.
: Good question, we assume it to be so, but that may be wrong. We moved the discussion Reply towards this question and away from functional description.

b) Page 18 -The information about genes found to be differentially expressed among morphs in your prior work should also be in the introduction, as it is background work that explains why you
F1000Research your prior work should also be in the introduction, as it is background work that explains why you took this transcriptomic approach. This can also be used to explain why you focused in on particular qPCR genes.
We added a sentence in the intro about the published papers, that this transcriptome made Reply: available. In those papers we focused on genes with putative craniofacial effects, though the focus in this study was broader.
c) A discussion of domestication related differences vs. benthic/limnetic differences should be included. I think the data from head gene expression is very interesting (Figs 5,S2) and really speaks to this question. d) In general, the role of stochastic evolutionary processes, and not just selection (artificial and natural) should be noted. For example, if the AC charr were simply taken from a stock with a different mtDNA haplotype then these differences in the mtDNA genome might not be adaptive, just random. If the AC fish has much higher mtDNA expression might this be simply a domestication issue and not indicative of selection in SB as stated? Finally, you find that not all mitochondrial transcripts (which are transcribed as a polycistronic transcript) are found at similar levels (Table 1) -what does this tell you about differential degradation/post-transcriptional processes? e) There is no discussion about the "Analyses of polymorphism in Arctic charr transcriptome" (Table 3,4,5), except for the mtDNA.
(c,d,e) Excellent suggestions. We added in the final discussion section few sentences on Reply: domesticated charr vs Benthic/limnetic. Unfortunately we do not have quantitative data on the phenotypes (head shape, and jaw) of the AC charr and acknowledge that we categorize it as limnetic based on general features.
We gladly added a sentence citing neutral forces, and are acutely aware that much of the divergence is likely due to history, drift etc. The domestication can certainly be the driver for the higher expression in AC -but we need transcriptomes from more populations/morphs to address that point. And yes, the variance in RNA levels from different parts of the mtDNA do indeed suggest differential half life of the various RNA species. Some are certainly degraded and others most probably actively utilized / protected. We decided not to follow that thought further though, as the MS already consists of quite a few threads already.
We also added sentences on the genetic polymorphism, before focusing more on the mtDNA. The main reason we dont want to elaborate to much on the SNPs is that we feel these data are mainly for generating hypotheses, and that more work is needed to substantiate SNPs and study their distribution in other populations.

Minor Comments
Introduction: a) Paragraph 3 -"Furthermore, recent estimates from the rainbow trout….by utilizing multiple data sources the genome assembly problem of this family can be solved". I am not sure how this statement is relevant to this particular study. This and the following statement seem more appropriate for the methods/discussion to me.
F1000Research appropriate for the methods/discussion to me.
We deleted this sentence and simplified the paragraph.

Reply:
b) The morphs being discussed should be clarified throughout the paper. For example, the authors often state "among morphs/among charr populations" but it is not clear which of the many morphs they are referring to (e.g. Paragraph 5, first sentence on allozymes and mtDNA and later sentence on MCHIIa -do you mean all 4 morphs of specific 2-way comparisons? Are some morphs more differentiated than others?) We tried to clarify this in various places in the manuscript, but in some cases we refer to Reply: morphs in general. Genetic separation can be estimated with Fst values either between pairs or over a larger set of groups (populations, morphs). In the intro we cite the work done to date in Iceland, which highlights the need for more pop. genetic analyses.

Methods:
a) The authors should note why they did not use the PI (large piscivorous) morph in any qPCR studies (in the methods or discussion) as this would be a nice morph to use in their tests for parallelism.
The PI charr is very rare in the lake and hard to catch. We later captured few sexually Reply: mature individuals, and generated couple of families, that were used for one study (Ahi et al Evodevo 2015). b) Page 5 (last paragraph) -the methods used to remove particular variants needs to be clarified. In particular, why the assumptions used to remove variants are valid by referencing past studies.
Many of the principles are common to most pipelines for removing spurious variants. In Reply: addition we applied filters necessitated by the properties of our dataset (pool of individuals), the mapping to an outgroup and paralogs due to salmonid genome complexity.  Figure 7 -It is not clear to me which variant is present in which morph. Adding the nucleotide to the x-axis (i.e. frequency of m1829G for B) would make this figure easier to quickly interpret. The "A.charr_WT" and "A.charr_M" should also be defined in the legend and it would be more appropriate to use scientific names for all species. Now fixed Reply:

Discussion:
F1000Research a) Discussion of reference 32 -The discussion of reference 32 is not put into the proper context. Figure 6 of this paper (Berthelot et al. 2014) shows that there are many genes that have no correlation among expression patterns and/or differences in expression levels (1573, 1248, and 1895=4716 paralog pairs), and that together these represent more than the 1,407 correlated/similar expression level paralogs. This section of the discussion needs to be modified.
Really valuable point, that we are especially grateful for. That we have added this fact to the Reply: intro and altered our interpretations in the discussion.
b) The Norman et al. (2014) paper should be mentioned earlier -if this is available why was it not used for their analyses? As well, the last sentence in this paragraph can be cut as it is evident.
The Norman papers are now presented more clearly in the intro. There are historical Reply: reasons for not including their data in our analyses, we had completed the analyses for this manuscript when they became available and have since then focused our data analyses efforts on another transcriptome generated in the lab (with longer reads).

Salvelinus alpinus
The work is founded on the solid premise that rapidly evolving phenotypes in nature can be underpinned by changes at the transcriptome level. The model system here is Arctic charr populations that have evolved (since the last ice age) major differences in phenotypes along the 'benthic' -'limnetic' axis, with strong differences in head morphology linked to feeding specializations. The work provides an extensive analysis of transcriptome and genetic differences between different morphs and populations. It is interesting, generally well-written and has merit on many levels. It is also rather hard going, since so much ground is covered on diverse areas. The study also comes with a large number of caveats, of which the authors are undoubtedly aware. Overall though, I am supportive of this work, as it represents one of the most detailed analyses of molecular mechanisms linked to rapid phenotypic evolution in Arctic charr. I see it as a great start point for future work and a source of several new findings and hypotheses. I suggest that F1000Research 1.

2.
it as a great start point for future work and a source of several new findings and hypotheses. I suggest that the paper be indexed in F1000 Research as long as its caveats are transparent and the authors address my comments.
I list below a number of suggestions that may help the authors improve the work, or that at least highlight study limitations for the benefit of interested readers. I also provide a number of minor comments and suggestions, which should help improve the manuscript more incrementally.

Main comments & caveats
RNAseq study design. I sympathize with the fact that the authors are trying to publish Illumina data that was generated in 2009, since (obviously) the technology has moved on greatly in the last 6 years, while its costs have been reduced dramatically. Adding to this is the fact that the authors are using a particularly complex transcriptome in terms of high content of similar paralogues (and expressed transposable elements), without a reference sequence for mapping in their species. I accept the author's argument that it is more sensible to map against a closely related species with the sequence data rather than to try and create a assembly from 36bp reads. I also de novo believe it is sensible to pool read counts for putative paralogous contigs in this study, since the short read length ablates any ability to separate paralogous differences in expression (yet does not preclude the generation of useful hypotheses about putative gene expression differences among morphs).
However, I do question whether the use of Atlantic salmon EST contigs is the best approach here. Firstly, reference assemblies for both Atlantic salmon and rainbow trout are now available, which distinguish paralogous variation. More importantly, using these reference genome data would provide certainty that reads are being mapped to exons from single genes, whereas many of the ESTs will provide a fragmented representation of exon sequences, presumably relying on annotation to piece them back into 'genes' . In addition, paired 100bp Ilumina reads are post hoc available at high coverage for Arctic charr (e.g. Norman . 2014), which could also be used to et al generate a specific reference transcriptome to map against in this study, although this might be underrepresented in terms of developmental genes as it is a gill study. Overall, I do wonder how much more information might have been gleaned from this dataset with a different mapping strategy?
With all the above said, I understand that the authors have built up a large study based around the original mapping to the salmon ESTs and that it would not be routine for them to repeat the study using better reference data. Furthermore, the approach used has definitely led to the generation of several valid hypotheses concerning the nature of gene expression and genetic differences among charr morphs, which have been followed up using independent approaches.
Methods "Biological Replication in RNAseq" -a general comment: obviously the design of the study is not optimal because biological variation within developmental stages is not considered in the statistics. Thus, the approach lacks power to detect differences when morph variation is restricted to different developmental stages. I wanted to explain my opinion (for the record) that the study design is nonetheless useful for identifying constitutive differences between morphs. This is especially true because gene expression variability is likely to be relatively low in embryonic stages (compared to a similar study design in adults at least). Further, the pooling of individuals will have helped to at least recapture some biological variation at different stages. Thus, as mentioned above, I see the author's use of RNAseq as a hypothesis-generating approach, which has been Introduction: "Examples of such a species complex are the finches of the Galapagos islands, . Grammatically cichlids in the African great lakes are exciting multi-species systems in this respect" -reads better: "Examples of such species complexes are provided by finches of the Galapagos islands, while cichlids of the African great lakes also provide an exciting multi-species system in the same respect" Introduction: "Some northern freshwater fish species exhibit frequent parallelism in trophic " change structures and life history and in several cases are they found as distinct resource morphs to " " …. are found as distinct resource morphs Introduction: " " change to "… in the development of ecological differences in tropic morphology ". trophic morphology ". This The family is estimated to be between 63.2 and 58.1 million years old information is not correct -it is correct to state that the age of the salmonid crown (based on the cited paper; different estimates exist in the literature, e.g. Macqueen and Johnston, 2014; ) is estimated at 63.2 and 58.1 million years old, but the family dates back Campbell . 2013 et al much further -to the origin of the WGD event in fact, which occurred more like 88-103 Ma (Macqueen and Johnston, 2014;Berthelot . 2014). Thus, the last common ancestor to extant et al salmonid species is what the authors are actually referring to in this sentence.
Introduction: "Furthermore, for data with short reads, mapping to a related reference ". While this sentence is genome/transcriptome is recommended over de novo assembly technically correct in the context of the work cited, I feel it is being used slightly out of context. For a start, what comprises a 'short read' is undefined. 36bp is short, but it is possible to get a sold reference transcriptome using 2*100bp, assuming the appropriate diversity of transcripts is represented and suitable depth is attained.
Introduction: " " change to " " nuclear genes, reveled both subtle nuclear genes, revealed both subtle Minor comment -AC, PL, LB and SB were already defined in introduction.
Methods: " " changed to " Fishing in Lake Thingvallavatn was with permissions Fishing in Lake ".

Thingvallavatn was done with permissions
Methods: " " change to " of differently expressed genes, we preformed clustering analyses …we " performed clustering analyses Results: "The most drastic changes were seen in processes related to glycolysis (GO:0006096, " change to "…. ". FDR = 0.0009), were the expression of 19 out of 25 genes where the expression Figure 7. What does the charr_WT vs. charr_M signify in the alignment data?
Discussion "We are interested in how predictable evolution is a the molecular level and if there " certain principles influence the rewiring of developmental and regulatory systems during evolution consider changing to "We are interested in the predictability of evolution at the molecular level, especially whether there exist principles that influence the rewiring of developmental and

Main comments & caveats
RNAseq study design. I sympathize with the fact that the authors are trying to publish Illumina data that was generated in 2009, since (obviously) the technology has moved on greatly in the last 6 years, while its costs have been reduced dramatically. Adding to this is the fact that the authors are using a particularly complex transcriptome in terms of high content of similar paralogues (and expressed transposable elements), without a reference sequence for mapping in their species. I accept the author's argument that it is more sensible to map against a closely related species with the sequence data rather than to try and create a de novo assembly from 36bp reads. I also believe it is sensible to pool read counts for putative paralogous contigs in this study, since the short read length ablates any ability to separate paralogous differences in expression (yet does not preclude the generation of useful hypotheses about putative gene expression differences among morphs).
However, I do question whether the use of Atlantic salmon EST contigs is the best approach here. Firstly, reference assemblies for both Atlantic salmon and rainbow trout are now available, which distinguish paralogous variation. More importantly, using these reference genome data would provide certainty that reads are being mapped to exons from single genes, whereas many of the ESTs will provide a fragmented representation of exon sequences, presumably relying on annotation to piece them back into 'genes' post hoc . In addition, paired 100bp Ilumina reads are available at high coverage for Arctic charr (e.g. Norman et al. 2014), which could also be used to generate a specific reference transcriptome to map against in this study, although this might be underrepresented in terms of developmental genes as it is a gill study. Overall, I do wonder how much more information might have been gleaned from this dataset with a different mapping strategy?
With all the above said, I understand that the authors have built up a large study based around the original mapping to the salmon ESTs and that it would not be routine for them to repeat the study using better reference data. Furthermore, the approach used has definitely led to the generation of several valid hypotheses concerning the nature of gene expression and genetic differences among 32 F1000Research using better reference data. Furthermore, the approach used has definitely led to the generation of several valid hypotheses concerning the nature of gene expression and genetic differences among charr morphs, which have been followed up using independent approaches.
: We thank the reviewer for excellent diagnosis and suggestions. The paper describes the (in Reply our humble opinion) most sensible summary of the data, as the writing of the paper started 2 years ago. We did map on the O.mykiss cDNA collection also, got similar results, but opted for reporting on the salmon data to avoid further extending an already long manuscript. We are currently analyzing DE and SNPs on a new assembly (100 bp PE reads -48 samples -3 morphsdevelopment), and may include a remapping of this dataset in that.
Methods "Biological Replication in RNAseq" -a general comment: obviously the design of the study is not optimal because biological variation within developmental stages is not considered in the statistics. Thus, the approach lacks power to detect differences when morph variation is restricted to different developmental stages. I wanted to explain my opinion (for the record) that the study design is nonetheless useful for identifying constitutive differences between morphs. This is especially true because gene expression variability is likely to be relatively low in embryonic stages (compared to a similar study design in adults at least). Further, the pooling of individuals will have helped to at least recapture some biological variation at different stages. Thus, as mentioned above, I see the author's use of RNAseq as a hypothesis-generating approach, which has been quite fruitful in identifying putative differences between different morphs.
: We appreciate the reviewers careful analyses of the study and approach. We tried to Reply emphasize the "hypothesis-generation" aspect during the rewrite.
Methods "QPCR study design". The authors adhere to the MIQE guidelines, but do not always follow the best approaches. Most pertinently, the authors use the 2−∆∆Ct method (assuming PCR efficiency of 2.0) despite having gone to the effort of gaining and reporting efficiencies for each assay, which can be as low as 1.72 for some genes. The effect of failing to incorporate differences in efficiency are highly established and this is likely to have affected the author's results. The authors should consider incorporating the effect of differences in efficiency into their analyses. This is likely to have some impact on the study conclusions in my opinion.
: Great point. The qPCR primer efficiencies more than 1.90 can be easily assumed as 2 Reply because of the negligible effects. Since we used LinReg software for efficiencies not the traditional method, it takes into account the efficiencies for each test for a given primer pair and discard those have different and lower efficiencies. However, the Natterin-like paralogues were below the cut-off. The statistical analyses were done on deltaCt values, prior to transformation based on efficiencies used for visualization. We now report the graphs of their expression adjusting for the lower efficiency, and state in the results "Note however, the efficiency of the primers for the nattl genes ranged from 1.72 to 1.77, which suggests this data should be interpreted with caution." Methods "Polymorphisms in charr transcriptome". While this is not exactly my area of expertise, I struggled to understand the methods behind filtering paralogous variants from SNPs in the data. The authors state "As the SNP analysis was done on individual contigs, differences among paralogs appear in the data. However, since each sample is a pool of few individuals, it is very unlikely that we have the same frequency of true SNPs in the samples. This property was used to remove variants that are most likely due to expressed paralogs". Can the authors please try to re-explain this in even simpler terms to help me get it? I don't see how this description leads to a robust identification of paralogous variation. Is there an underlying assumption of equal expression among paralogues? If so, this is likely to be routinely invalidated.
F1000Research among paralogues? If so, this is likely to be routinely invalidated.
: We acknowledge this part is a hard read. We rewrote this part of the methods. Here is Reply another summary. Reads from regions that are very similar in paralogous genes can map to both of them. Because we consider also reads that map to many contigs, some of the candidate variants will reflect sequence differences between paralogs, not polymorphism in either paralog. Next we deploy the population genetic argument, since we are sequencing RNA from 6 chromosomes in each sample, then it is very unlikely that a TRUE SNP will be at the same frequency in all of the 8 samples. But variants -that are due to differences bwn paralogs -are likely to be similar in frequency because they are unaffected by the population sampling. This filter is designed to toss those out.
To emphasize the objective is not to find differences between paralogs, but rather to enrich for true SNPs. This method will toss out many sites separating paralogous genes (but not all because some paralogous genes are differentially expressed between morphs or time points).
Methods "Verification of candidate SNPs". While it is good that the authors have attempted to verify SNPs identified from their RNAseq data, I don't believe the data is particularly well incorporated in the results section. It needs to be stated up front the extent to which the SNPs predicted from the RNAseq were independently verified. Also, the methods for this section can be improved, especially "we conducted genomic comparisons of the Salmon genome, ESTs and short contigs from the preliminary assembly of the Arctic charr transcriptome". None of this information is elaborated on -what is the preliminary assembly of the Arctic charr transcriptome? Which version of the salmon genome was used and how? Moreover, it would be useful to actually explain in the methods that the genotyping was done on a small number of SB, PL and PI morphs, rather than relying on the reader to extract all the required information from Table S2. I guess overall, the way this section is incorporated into the manuscript needs some thought in terms of improving the reader's experience. I struggled after reading it several times and am still not sure I have all the information I need.
: We fixed the methods section to accommodate both reviewers which brought up similar Reply points. We highlight the sampling (8 individuals of 3 morphs), and extend the description of the genomic comparisons. We also extend the discussion of those results.
Results. "Analyses of those reads require an Arctic charr genome sequence or transcriptome assembly from longer and paired end reads." As mentioned already, the latter is available to generate an Arctic charr transcriptome assembly to map against.
: Unfortunately the great Norman 2014 data Reply et al. (http://www.ncbi.nlm.nih.gov/pubmed/24368751) came to our attention after we had done these analyses, and started working on our new data (see above). Thus we opted for not redoing the whole analyses for this manuscript, but focus on the verification -and of course working on a new assembly using longer reads. Figure 3 and 4. The authors found that around half the genes studied were not differentially expressed among morphs by qPCR. Obviously this is quite a large number, but on closer inspection, I noticed that Ndub6, Ubl5 and parp6 were not even differentially expressed according to RNAseq. Thus, I am confused at the selection of genes from the RNAseq analysis for verification by qPCR. The authors should explain this selection more transparently and provide clearer indices of the correlation between RNAseq and qPCR results and associated discussion.

F1000Research
Introduction: "Furthermore, for data with short reads, mapping to a related reference genome/transcriptome is recommended over de novo assembly". While this sentence is technically correct in the context of the work cited, I feel it is being used slightly out of context. For a start, what comprises a 'short read' is undefined. 36bp is short, but it is possible to get a sold reference transcriptome using 2*100bp, assuming the appropriate diversity of transcripts is represented and suitable depth is attained.
: Great point, we opted for keeping the point (at this place in the ms) but changing the Reply wording to: In this study we opted to map the reads (36 bp) to a related reference genome/transcriptome {Vijay2013a}, instead of conducting de novo assembly.
Introduction: "nuclear genes, reveled both subtle" change to "nuclear genes, revealed both subtle" : Thanks fixed. Reply Minor comment -AC, PL, LB and SB were already defined in introduction.
: Thanks, removed this. Reply Methods: "Fishing in Lake Thingvallavatn was with permissions" changed to "Fishing in Lake Thingvallavatn was done with permissions".
: Ammended. Reply Methods: "of differently expressed genes, we preformed clustering analyses" change to "…we performed clustering analyses" : Thanks, fixed. Reply Results: "The most drastic changes were seen in processes related to glycolysis (GO:0006096, FDR = 0.0009), were the expression of 19 out of 25 genes" change to "…. where the expression".
: Thanks, fixed. Reply Figure 7. What does the charr_WT vs. charr_M signify in the alignment data?
: Designates the two alleles, the legend now makes this explicit. Reply Discussion "We are interested in how predictable evolution is a the molecular level and if there certain principles influence the rewiring of developmental and regulatory systems during evolution" consider changing to "We are interested in the predictability of evolution at the molecular level, especially whether there exist principles that influence the rewiring of developmental and regulatory systems".
: Thanks, excellent suggestion, included Reply Discussion. "Recent rainbow trout data shows most paralogs from the latest whole genome duplication event retain the same expression pattern32 indicating that this scenario is probably uncommon; hence it is of considerable interest when two paralogs show distinct expression patterns". I do not agree that it is of considerable interest when two paralogs show distinct expression patterns -I could list tens of examples for salmonids.
: Good point, we have revisited this interpretation (see also point by rev. 1). Reply

Conclusions "
The results suggest genetic and expression changes in multiple systems relate to divergence among populations." Change to "… associated with divergence among populations."