1. Introduction {#sec1}
===============

The deep-sea (\>1000 m depth) covers about 70% of the Earth\'s surface, representing one of the last large unexplored areas on the planet. Only within the last few decades the technology has advanced sufficiently to reach the deep-sea effectively, revealing unexpected high levels of biodiversity and extremely diverse habitats (canyons, cold seeps, hydrothermal vents, deep-water coral reefs, mud volcanoes, seamounts, and trenches) of significant conservation interest and potential high economic values. Deep-sea environments are characterized by extremely high hydrostatic pressures (1 MPa every 100 m), lack of light, and low temperatures (down to 1-2°C). Therefore, fish as well as any other organism living in the deep-sea had to adapt to tolerate conditions of this extreme habitat \[[@B1]\].

First studies on adaptation to high pressure and low temperatures are dated back in the '70s and they report comparison of common proteins present in shallow and deep-water fishes \[[@B2], [@B3]\]. Key enzymes in muscle tissues that exhibit adaptive differences among species at different depths are the lactate dehydrogenase (LDH) and malate dehydrogenase (MDH) presenting differences in structural stability (reviews in \[[@B4]--[@B6]\]). More recent studies on evolutionary adaptation of functional genes to high pressure report unique amino acid substitutions in *α*-skeletal actin and myosin heavy chain (MyHC) proteins in deep-sea fishes \[[@B7]--[@B10]\]. For deep-sea species inhabiting hydrothermal vents and cold seeps, environments characterized by high pressure, chronic hypoxia, and high concentrations of toxic compounds, molecular and functional adaptation of hemoglobins (Hbs) are reviewed in Hourdez and Weber \[[@B11]\]. Despite these studies, our knowledge on wide scale gene expression patterns in deep-sea fish remains elusive.

The black scabbardfish,*Aphanopus carbo* (Lowe 1839), is a bathypelagic species belonging to the Trichiuridae family and is distributed in temperate-cold Atlantic waters at depths between 200 and 1800 m \[[@B12], [@B13]\].*A. carbo* represents a commercially valuable species for several regions of the Iberic peninsula, especially in Madeira where catches have reached up to 1000 tons per year \[[@B14]\] amounting to ca. 55% of the total landings. Recently this species has become increasingly targeted by Portuguese, French, and Irish fishing fleets (\[[@B15]\] and literature therein) and fishery data have shown a constant decline in population \[[@B16]\]. The information available on the biology, maturity, spawning, and growth of this species \[[@B17], [@B18]\] is scattered. Recent studies are reporting a panmictic distribution of this species in the NE Atlantic with multiple breeding sites at low latitudes \[[@B19]\]. It is also worth mentioning that, in southern locations, this species lives in sympatry with*A. intermedius*, a close related species with very similar morphology \[[@B20]\], therefore attracting interest for evolutionary studies.

High-throughput sequencing approaches applied to transcriptomics now provide a global perspective on taxonomic and functional profiling of genes expectedly expressed under the influence of environmental conditions in which these organisms live. Also known as next-generation sequencing, these techniques allow for a massive characterization of expressed sequence tags (ESTs) providing an overview of those genes expressed in a given tissue at any given time \[[@B22]\]. In silico analyses of massive gene libraries may serve several interests among others. For instances, from discovery and identification of new genes, characterization of gene expression, to development of novel genetic markers for quantitative trait locus (QTL), and population of genomic analyses. The breadth of next-generation sequencing applications extends over a variety of biological questions including those addressing pertinent questions regarding a species\' ecology, life history, and evolution \[[@B21], [@B23]\].

Previous studies regarding transcriptome sequencing and gene expression studies in deep-water species were mostly limited to hydrothermal vents invertebrates \[[@B24]\], microbial communities in hydrothermal plumes \[[@B25]\], deep-sediments \[[@B26]\], and in the water column \[[@B27]\] leaving vertebrates species virtually under-represented. The present work represents a pioneer study for deep-sea fishes providing new insights into the role of differential gene expression on the environmental adaptation of deep-sea black scabbardfish.

Here we describe the assembly and annotation of the transcriptome of*A. carbo* obtained by sequencing mRNA libraries of six tissues (spleen, brain, heart, gonads, liver, and muscle) and explore functional genes whose sequence might be associated to depth adaptation. Additionally, we tested the correlation of selected candidate genes comparing the number of contigs against the gene expression normalized to a relative value of 1.0, as determined by qPCR amplification. Furthermore, the screening of the libraries allowed the identification of new EST-simple-sequence repeats (SSRs) and the design of primer pairs suitable for population genetic studies as well as tagging and mapping of genes.

2. Methods {#sec2}
==========

2.1. Fish and Tissue Samples {#sec2.1}
----------------------------

Specimens of*Aphanopus carbo* (two males and two females) were collected in 2009 onboard of the RV "*Arquipelago.*"*A. carbo* were fished at depth range 1100--1250 m using deep-water long-lines in proximity of the Condor Seamount, located approximately at 15 nm SW of the island of Fayal (Azores, Portugal). The four specimens used in this study were caught on the same longline set and once onboard, the freshly caught animals were dissected and portions (or complete organs) of spleen, brain, heart, gonads, liver, and muscle tissues were preserved in formamide solution and kept at −20°C until RNA extraction was performed.

To validate the correct identification of the species, a small portion of muscle tissue was also preserved in 95% ethanol for molecular screening following Stefanni et al. \[[@B28]\] protocols.

2.2. RNA Extraction and Sequencing {#sec2.2}
----------------------------------

Total RNA was extracted from 20 to 40 mg of each of the six preserved tissues of a pool of four*A. carbo* individualsusing the RiboPure kit (Ambion, Applied Biosystems). Quantity and purity of the RNA was determined on a 1.4% agarose-MOPS-formaldehyde denaturing gels and by assessing the *A* ~260/280~ and *A* ~260/230~ ratios using the NanoVue spectrophotometer (GE Healthcare). Poly-A RNA was extracted from 15 *μ*g of each total RNA sample using the Poly(A)Purist mRNA Purification kit according to manufacturer\'s instructions (Ambion, Applied Biosystems). mRNAs were transcribed into cDNA utilizing Mint-2 cDNA synthesis kit (Evrogen, <http://www.evrogen.com/>) according to manufacturer\'s instructions for NGS platforms. Six cDNA libraries were constructed from mRNA of individual pools of tissues and sequenced in a single 454 GS FLX Titanium run. Each of the cDNA libraries was characterised by unique sequence tags (MIDs) that allowed to trace back the sequences generated from single tissues after assembly.

cDNAs were sheared by nebulization to yield random fragments approximately 500--800 bp in length, by applying 30 psi (2,1 bar) of nitrogen for 1 minute on 4 *μ*g of each library. The distribution of fragments was verified on a BioAnalyser DNA 7500 LabChip (Agilent Technologies). The fragmented cDNA sample was end-repaired with T4 DNA polymerase and T4 polynucleotide kinase and adaptor sequences ligated according to the manufacturer\'s instructions \[[@B29]\]. The fragments were immobilized onto streptavidin beads and nick-repaired with*Bst* polymerase. The cDNA fragments were denaturated with alkali to yield single stranded cDNA (sscDNA) library. Quality of the library was assayed on a BioAnalyser RNA 6000 Pico LabChip (Agilent Technologies) and quantity measured by spectrofluorimetry with the Quant-iT RiboGreen RNA Assay kit (Invitrogen). A titration was set up at 1, 2, 4, and 8 copies per bead (cpb) in the clonal amplification by emulsion PCR to optimize yield and sequence quality. The percent enrichment of beads carrying the sscDNA was determined and the amount of library input calculated to 18%. A large scale emulsion PCR was set-up based on the previous value and sequenced at Biocant (Cantanhede, Portugal) using the 454 GS FLX Titanium pyrosequencing on a full 70X75 PicoTiterPlate, according to manufacturer\'s instructions (Roche).

2.3. Bioinformatic Analyses {#sec2.3}
---------------------------

High quality reads were assembled using Newbler ver. 2.6 (Roche 454) sequence analysis software. All reads were identified and grouped by their unique MIDs to the tissue of origin. Trimming and masking the polyAs was a common procedure for the assembling tool.

The assemblage is characterized by read overlaps and multiple alignments made in nucleotide space. Consensus base-calling and quality value determination for contigs are performed in flow space. The use of flow space in determining the properties of the consensus sequence results in an improved accuracy for the final base-calls. The implementation of this software was performed using default parameters. Assembled contigs were annotated through sequence similarity searches against the National Centre for Biotechnology Information (NCBI) nonredundant (nr) protein database using the BLASTx \[[@B30]\] with a cut-off criterion of an expect-value (*e*-value) \< 10^−6^. The contigs that did not find a hit were further processed with ESTScan (<http://www.ch.embnet.org/software/ESTScan2.html>). The two assemblages of amino acid sequences, resulting from the BLASTx searches at high level of stringency and the ESTScan, were processed by InterProScan for functional annotation of transcripts applying the function for the mapping of gene ontology (GO) terms. The GO method classifies genes within a hierarchy using a systematic nomenclature of attributes that can be assigned to all gene products independently from the organism of origin. To reduce the redundancy in the consensus sequences which correspond to the same gene we used BLASTClust to detect similar assemblies with 95% identity and 90% coverage. All the results from both assemblage methods were loaded into a SQL database developed for this purpose.

To validate the accuracy of the assembly, the resulting contigs were compared to previously sequenced transcriptomes of 6 teleosts including*Danio rerio*,*Gasterosteus aculeatus*,*Oreochromis niloticus*,*Oryzias latipes*,*Takifugu rubripes*, and*Tetraodon nigroviridis*, using tBLASTn \[[@B30]\] to find protein homologs at two levels of stringency (*e*-value \< 10^−3^ and *e*-value \< 10^−10^).

To identify protein conserved domain specific for each tissue analysed a new annotation was performed with Hmmer against the Pfam database (ver. 25.0). Protein domain representativeness for each tissue was obtained comparing protein domain abundance in a particular tissue versus all the tissues compiled together using a hypergeometric test.

2.4. cDNA Synthesis and qPCR Validation Tests {#sec2.4}
---------------------------------------------

Fresh cDNA was synthesised from the six mRNAs that were used for pyrosequencing, cDNA synthesis was performed using primers with oligo(dT) and the ThermoScript RT-PCR System (Invitrogen) following the manufacturer\'s instructions.

A set of 28 genes were selected including candidates that were tissue-specific and genes that were encountered in the tissue expressed at similar as well as at different amounts in all the six libraries, with the aim of covering most of the possible expression scenarios within the dataset. Frequencies of contigs for all candidates genes in the mRNA libraries were obtained by detecting orthologous gene sequences using the BLAST tool included in the*A. carbo* database.

For the design of all qPCR primers ([Table 1](#tab1){ref-type="table"}) we used the web interface NCBI Primer-BLAST (<http://blast.ncbi.nlm.nih.gov/>). Alignments of the sequences provided by the output from the internal blast search were used to select all primer sets.

Gene expression, calculated as relative expression, was determined by means of real-time PCR using the CFX96 (Bio-Rad). Primer concentrations and sample dilutions were optimized to meet highest efficiency in the PCR reaction in a total reaction volume of 20 *μ*L. Fluorescent signal was determined by the addition of SsoFast EvaGreen Supermix (Bio-Rad) which was included in the cocktail accordingly to manufacturer\'s instructions. Baseline and threshold cycle were always set to automatic in the sequence detection software, CFX Manager (Bio-Rad). All plates contained a "no template control" (NTC) and each sample was tested in duplicate. Cycling conditions for gene amplifications were 95°C for 3 min followed by 35 cycles of 95°C for 10 s, 56°C for 10 s, and 68°C for 15 s. An additional protocol for melting curves analysis included a cycle at 95°C followed by a progressive reading of fluorescence for every cycle from 65°C to 95°C for 5 s at intervals of 0.5°C. Gene expression, normalized to a relative value of 1.0 for all the genes selected and for each tissue, was compared to the contigs frequencies generated by the assembly platform to determine the significance of correlation between qPCRs values and 454 sequencing reads from unnormalized cDNA libraries.

2.5. Characterization of Depth-Related Functional Genes {#sec2.5}
-------------------------------------------------------

The predicted amino acid sequences of functional genes of*A. carbo*possibly related to depth adaptations were compared to those deposited in NCBI database searching for homologies. We aligned and compared translated sequences of lactate dehydrogenase (LDH-A and LDH-B), cytosolic malate dehydrogenase (MDHc), hemoglobins (Hb-A and Hb-B), actin (ACTA1), and myosin heavy chain (MyHC). Protein and nucleotide sequences were aligned using Clustal X \[[@B31]\] while sequence analysis and phylogenetic inferences were performed using CLC Main Workbench (v. 6.8.2, CLC Bio). The neighbor-joining (NJ) algorithm \[[@B32]\] was implemented to construct a phylogenetic tree using HKY substitution model and attributing a gap penalty of 10. The support for internal branches was assessed using the bootstrap \[[@B33]\] with 1000 replicates.

Nucleotide alignments and ML trees built, implementing the most appropriate substitution model under the Akaike Information Criterion (AIC), were used in the program "codeml" in PAML 4 \[[@B34]\] to assess selective pressure on those genes for which the complete sequences was available. Positive (or negative) selected sites were defined by the ratio between nonsynonymous*versus* synonymous substitutions (*dN*/*dS* or *ω*). Two models were tested comparatively: M1 which groups codons in two classes (*ω* \< 1 and *ω* = 1), clustering sites under negative or neutral selection; and M2 which groups codons in three classes (*ω* \< 1, *ω* = 1 and *ω* \> 1), adding a cluster for sites under positive selection to the ones defined in M1. Probabilistic measures of how well these models fits the evolutionary relationship of individual genes were calculated from the likelihood values of fitted models and the number of "free parameters" for all genes.

Protein stability was estimated using a virtual quantification software \[[@B35]\] calculating Gibbs free energy in terms of kinetic and thermodynamic quantities taking into account each amino acid contribution for the maintenance of the native structure of the protein. For a protein to maintain its stability there is a need of sufficient hydrophobic residues which will utilize free energy to guide proper folding \[[@B35]\].

2.6. EST-SSR Resources for Population Genetics {#sec2.6}
----------------------------------------------

Among various molecular markers, simple sequence repeats (SSRs) are highly polymorphic, easier to develop, and very useful for researches such as genetic diversity assessment. Therefore*A. carbo* database was further used to detect such regions in the transcriptome sequencing data and provide a list of combination of primer sets on flanking regions that could be used for population genetic studies. To identify EST-SSRs, all the contigs were searched using MISA \[[@B36]\] and for primer design we used Primer3 \[[@B37]\]. The algorithm of the SSR finder identifies a good quality repeat when one locus is present with adjacent loci at an up or downstream distance higher than 100 bp and parameters were set to locate a minimum of 20 bp sequence repeats: di-mers (x12), 3-mers (x8), 4-mers (x5), 5-mers (x5), and 6-mers (x5). Primer design was performed setting parameters of a minimum size of 20 bp and melting temperatures of 60°C.

3. Results and Discussion {#sec3}
=========================

3.1. Sequences Assemblage and Functional Annotations {#sec3.1}
----------------------------------------------------

After sequence trimming, a total of 544,491 high quality reads were produced with an average length of 237 bp corresponding to 129.5 Mb. A total of 8,319 contigs were assembled with Newbler (ver 2.6) as high quality consensus sequences without the presence of singletons. A summary of EST data for each of the six tissues is reported in [Table 2](#tab2){ref-type="table"}. A total of 2,440 assembled contigs were annotated against the NCBI nonredundant protein database at the cut-off (*e*-value) \< 10^−6^. Additional 3,843 contigs with no protein matches were further processed with ESTScan to find 2,715 homologous proteins. All the contigs were then processed by InterProScan for functional annotation of transcripts and mapping functional information to gene ontology (GO) terms. Of the 5,155 amino-acid sequences (out of the 6,283 totally sequenced), 1,728 could be annotated within the GO hierarchy ([Figure 1](#fig1){ref-type="fig"}), while 2,395 could be annotated accordingly to InterProScan. The complete GO mapping of*A. carbo* for individual tissues can be accessed on the dedicated online database (<http://transcriptomics.biocant.pt/AphanopusCarbo/>). The largest proportions of GO functional categories are of similar proportion to the unnormalized libraries of Tilapia \[[@B38]\]. In the Cellular Components GO group ([Figure 1(a)](#fig1){ref-type="fig"}), the genes involved in cell and cell part functional categories corresponded to the largest percentage of the pie chart (28.15% each). In the Molecular Function GO group ([Figure 1(b)](#fig1){ref-type="fig"}), almost half of the total genes were involved in binding (47.40%), followed by those related to catalytic activity (27.53%). Finally, for the Biological Processes GO group ([Figure 1(c)](#fig1){ref-type="fig"}), the largest portion of the genes were involved in metabolic processes (33.45%), tightly followed by the category of genes linked to the Cellular Process (31.93%).

To further validate the accuracy of the A. carbo assembly, we compared our dataset to the transcriptomes of six other prior sequenced teleosts including*Danio rerio*,*Gasterosteus aculeatus*,*Oreochromis niloticus*,*Oryzias latipes*,*Takifugu rubripes*, and*Tetraodon nigroviridis*. Similarity searches were performed comparing our assembled contigs against each of the available transcriptomes using tBLASTn \[[@B30]\]. The results highlighted strong similarity between the transcriptomes of*A. carbo* and other teleosts indicating that about 30% of the total contigs matched protein homologs (*e*-values \<10^−3^ ranging between 2,268 and 2,457) and that our assembly presents the highest similarity (by little margin) with the transcriptome of Gasterosteus aculeatus (*e*-value \< 10^−10^ = 2,223) ([Table 3](#tab3){ref-type="table"}).

A total of 1639 (20%) transcripts were annotated by Pfam protein domains matches and the set of protein domains characterizing each single tissue was identified by hypergeometric test ([Table 4](#tab4){ref-type="table"}).

3.2. qPCR Assays and Validation Tests {#sec3.2}
-------------------------------------

Quantitative PCR assays were performed using cDNA samples at 1 : 400 dilution (up to 100 ng approximately in the qPCR master mix), conditions at which the PCR resulted to be more efficient. We tested seven series of dilutions ranging from 1 : 50 to 1 : 1000. Based on optimized qPCR conditions; all targeted cDNAs were tested across all tissues for the complete set of genes selection used in this study. Normalized expression to a relative value of 1.0 for 28 genes is shown in [Figure 2](#fig2){ref-type="fig"} (see: Figure S1 in Supplementary Material available online at <http://dx.doi.org/10.1155/2014/267482>).

Statistical tests were performed to identify the correlation between the mean value of normalized expression and a relative value of 1.0, as determined by qPCR against contig frequencies using Pearson\'s *r* coefficient and its significance ([Figure 3](#fig3){ref-type="fig"}). Most of the genes showed highly significant correlation; however one case indicated a reduced significance (elongation factor) and a few other cases resulted to be not significant (e.g., Ras-related GTP-binding protein, Basic Transcription Factor 3, Cu/Zn Superoxide Dismutase and 2-Cys Peroxiredoxin). Reduced or lack of significance appeared to be related to a low number of contigs from the 454 sequencing ([Table 5](#tab5){ref-type="table"}), which might be regarded as an indication of failure to reach library saturation. Low values in contig frequency should correspond to a reduced level of expression as the libraries are nonnormalized; therefore caution should be used when exploring 454 outputs below a certain threshold. Our data derived from the reading of a single plate were insufficient for systematical comparison of all the genes, although the majority of the genes show a high correlation between the 454 data and the qPCR. However, it should be emphasized that this exercise was not meant to be a surrogate to the qPCR approach for the study of gene expression but more a proxy for a preliminary screening of differentially expressed genes in multiple libraries.

3.3. Candidate Genes Associated to Depth {#sec3.3}
----------------------------------------

The predicted amino acid sequences of*A. carbo* functional genes putatively related to depth adaptations were compared to those deposited in NCBI database for homologies searches. Complete alignments of orthologous protein sequences from the set of functional genes that have been reported as responsive to depth adaptations included representatives of Teleosts belonging to several families ([Table 6](#tab6){ref-type="table"}). Exploring the levels of similarities expressed as percentage of amino acid identity or number of amino acid that differ among sequences between*A. carbo* and other fish species ([Table 7](#tab7){ref-type="table"}) brings evidence supporting the fact that deep-living as well as polar species are not the most similar to the black scabbard fish. We attempted to reconstruct the phylogenies for those species using amino acid sequences of functional genes (LDH-A, LDH-B, MDHc-A, Hb-A and Hb-B) and corresponding mtCOI nucleotide sequences ([Figure 4](#fig4){ref-type="fig"}). The resulting trees showed similar topologies, suggesting that the signals embedded in these functional genes reflect evolutionary divergence among taxa rather than any enzyme adaptations relationships.

The lactate dehydrogenase A (LDH-A: isotig01401) was encountered abundantly in the muscle tissue ([Table 5](#tab5){ref-type="table"}, [Figure 2](#fig2){ref-type="fig"}) and its comparison with other orthologous sequences differs for a number of amino acids between 16 and 44. One indel at position 75 had to be added to the two polar species to obtain a complete alignment ([Table 7](#tab7){ref-type="table"}). The percentage of identity was higher with the strictly marine species ranging between 95.2 and 93.7%. Unfortunately, no sequences from neither deep-sea nor abyssal fishes were available for comparison. On the other hand, lactate dehydrogenase B (LDH-B: isotig01423) in*A. carbo* was highly represented in the heart while very scarcely recorded in the other tissues (Table S1). The largest sequence divergence in terms of amino acid substitutions ranges between 47 and 55 when compared with deep-water macrourids, and similarly to LDH-A, an indel had to be added in the sequences of the gadiform species at position 76 to obtain a complete alignment ([Table 7](#tab7){ref-type="table"}). Previous studies on benthopelagic fish report a significant decline of LDH activities with increasing water depth \[[@B39]\]. However, it is known that variation in enzyme activities of fish at a given depth is influenced by feeding behaviour and locomotory modes \[[@B40], [@B41]\].

Cytosolic malate dehygrogenase A (MDHc-A: isotig01010) was detected in most of the tissues, with the exception of the muscle (Table S1). An isoform of MDHc was also detected but its sequence covered only the last 110 amino acids (isotig01731). However the type of substitutions and tissue expression pattern (Table S1) according to Merrit and Quattro \[[@B42]\], lead us to believe that this cytosolic isoform is MDHc-B. It has been reported that the two teleost MDHc isozymes are the products of a gene duplication event after the separation of teleosts and tetrapods (see \[[@B42]\] and references therein) although the exact timing of this duplication has not been inferred.

Further comparison analyses were carried out only for MDHc-A, with orthologous sequences obtained from NCBI database where only shallow-water fish sequences were available. The complete alignment did not include any indel, and sequences differed from*A. carbo* between 15 and 35 amino acid substitutions, representing a percentage of identity ranging from 95.5 to 89.5%.

Two *α*-skeletal actins were detected from the*A. carbo* database: the isoform 1 (isotig01459), expressed in shallow-water species and probably not functioning in abyssal fishes, and the isoform 2a (isotig01561). The mutations specific for*A. carbo* in the actin 2a were the substitutions Asp3Glu, Ala155Ser (a mutation also common with the isoform 2b \[[@B7]\]), Ser234Val, Val165Ile, Leu261Val, and finally Ala278Thr. This type of isoform was reported in other deep-sea species as*Coryphaenoides acrolepis*,*C. cinereus*,*C. yaquinae,* and*C. armarus* \[[@B7]\]. The isoform 2b, so far reported only in abyssal species, was not detected in*A. carbo*.

Actins 1 and 2a were significantly expressed in muscle tissue with evidence of its presence also in the heart (Table S1). Direct comparison of*A. carbo* actin 1 with orthologous sequences reported for other marine and deep-water species differed by just 1 amino acid at position 3 ([Table 7](#tab7){ref-type="table"}). The number of substitutions increases to 2 when this sequence is compared to the isoform 2a of*Coryphaenoides*species (substitution at position 155) and to 5 amino acids replacement ([Table 7](#tab7){ref-type="table"}) when compared to the isoform 2b, two of which are unique (positions 116 and 137) \[[@B7]\]. Actin has a function in polymerization of G-actin to F-actin in neutral salts. While the volume is increased following polymerization, the reaction is strongly affected by high pressure \[[@B43]\]. It is reported that the substitutions of Val54Ala or Leu67Pro reduced the volume change associated with actin polymerization \[[@B7]\].

The assemblage of the myosin heavy chain protein (MyHC: isotig02568 and isotig1394) obtained in*A. carbo* was incomplete (All_gs454_000756 and All_gs454_000001) containing a gap of 711 amino acids at the positions from 171 to 882 of the 1933 AA of the complete sequence. Unfortunately, within this gap there are the two loop regions with characteristic structures that are uniquely reported for deep-sea fish: loop-2 region is shorter and the loop-1 region has a proresidue \[[@B9]\]. MyHC was almost uniquely found in muscle tissue (Table S1) and for comparison analyses with orthologous genes from other fish we only used the last 305 AA (All_gs454_000184) of the complete sequence. The lowest number of amino acid substitutions ranged between 53 and 69 (corresponding to sequence identity between 94.95 to 93.43%) when the sequence from*A. carbo*is compared to its relative shallow water marine and freshwater species, while this value increases up to 92 amino acid substitutions (and corresponding sequence identity of 91.25%) when compared to its relatives from deep water or polar regions ([Table 7](#tab7){ref-type="table"}).

Variation in globin sequences was analysed exploring the *α* ~2~- and *β* ~2~-chains (Hb-A: isotig00247 and Hb-B: isotig00163), whose relative expressions were detected more abundantly in the spleen, followed by liver, heart, brain, and muscle and virtually absent in the gonads ([Figure 2](#fig2){ref-type="fig"}, [Table 5](#tab5){ref-type="table"}). Comparison analyses with orthologous sequences of deep-sea gadiforms and a notothenoid indicate that the number of amino acid substitutions ranged between 39 and 53 (Hb-A) and between 31 and 42 (Hb-B). The complete alignments included the additional single indel for the *α* ~2~-chain of*Notothenia angustata* at position 102. Notothenioids acquired a completely different globin genotype with respect to other teleostean groups. The Antarctic ichthyofauna (dominated by a single taxonomically uniform group) lost its globin multiplicity in correlation with temperature stability. On the other hand, for the Arctic ichthyofauna it may have been advantageous to maintain a multiple globin system, helping to deal with environmental changes and metabolic demands \[[@B44]\].

Selective pressure in the site-by-site patterns among species was evaluated for LDH-A, -B, MDHc, and ACTA1 (isoform 1). There was no evidence of positive selection at the nucleotide site level in any of those genes whose global *ω* value was very low in all cases (from 0 in ACTA1 and 0.032 in LDH-B) (Table S2). The proportion of sites supporting positive selection, the model M2, was null (LDH-A and -B) or extremely low (0.3% in MDHc and 1.6% in ACTA1) (Table S2). These results indicate that the evolution of these four genes in teleosts is constrained by very stringent selective pressure (model probabilities for M1*versus* M2 ranging from 87% in MDHc and 88% in all others).

In terms of protein stability calculated as Gibbs free energy and taking into consideration kinetic and thermodynamic quantities, we explored only the functional genes for which the complete sequences were obtained ([Figure 5](#fig5){ref-type="fig"}). In this context, protein stability is defined by the ability of a protein to retain its structural conformation or its activity when subjected to physical or chemical manipulations; therefore the energy consumed by activation to promoted folding has to be compensated by thermal stability provided by the energy of denaturation \[[@B35]\]. Hence, protein stability is quantitatively calculated by the standard Gibbs energy change (Δ*G*), allowing comparison of stabilities for different proteins \[[@B45]\].

The normalised difference in protein stability (Δ*G*) of lactate dehydrogenase A (LDH-A) was lower in*A. carbo* compared to fish from polar waters and tropical marine ([Figure 5](#fig5){ref-type="fig"}). This gradual variation was primarily dictated by the lower contributions of leucine and lysine in the hydrophobic trend of kinetic calculation in the*A. carbo*sequence (Table S3). The substantial drop in Δ*G* of LDH-B in*A. carbo* compared to the orthologous sequences of polar, abyssal, and shallower marine fishes ([Figure 5](#fig5){ref-type="fig"}) was due to a larger contribution of glutamic acid (Table S3) promoting the thermal denaturation of the protein. Other studies investigating kinetic, physical properties and ability to withstand high pressure of LDH-B of two gadiformes support our calculations showing that the enzyme from the deep-sea species has a significant increased tolerance to pressure and higher thermostability \[[@B6]\]. To provide protection to*A. carbo* LDH-B from pressure and temperature denaturation, osmolytes might play an essential role. Experiments adding Trimethylamine-N-oxide (TMAO) to samples resulted in substantial increment of activity of LDH-B under conditions at which the enzyme was previously sensitive \[[@B6]\].

In the cytosolic malate dehygrogenase A (MDHc-A) there was a sharp decrease in protein stability of marine fishes going from shallower to deeper waters (no sequences of abyssal species were available). Lowest values are encountered in the sequences of the two euryhaline species included in the comparison ([Figure 5](#fig5){ref-type="fig"}). This progressive decrement in stability was associated to the larger contribution of aspartic acid (Table S3), also a promoter of thermal denaturation of the protein. The responses to pressure and temperature of soluble enzymes like LDH and MDH differ adaptively among species found at different depths.

The isoforms 1 and 2a of *α*-skeletal actin of*A. carbo* showed very similar numerical values of Δ*G* compared to their homologues of deep-water (actin 1) and deep and abyssal water (actin 2a). Protein stability drops substantially when these two isoforms were compared with actin 2b, only present in abyssal species ([Figure 5](#fig5){ref-type="fig"}). Such reduction in stability was due to a higher number of hydrophobic amino acids contributing to the thermodynamic calculation (Table S3). Morita \[[@B7]\] reported that the substitutions of Gln137Lys and Ala155Ser generate a mechanism for stabilizing enzyme-substrate interactions under high pressure.

Comparing the stability of globin, *α* ~2~-chain (Hb-A), in shallow, deep-water, and polar fish, the larger Δ*G* was found in polar and deep-water cods, followed by*A. carbo* before dropping to negative values for the shallower representatives of Antarctic and temperate environments ([Figure 5](#fig5){ref-type="fig"}). Stability of protein in*A. carbo* and its deep-water relatives was linked to larger contributions by leucine in the kinetic calculations promoting folding and thermal stability (Table S3). The stability as Δ*G* of globin sequences of the *β* ~2~-chains (Hb-B) showed a very similar trend as in Hb-A (although there were no negative values), and the two deep-living species had very similar protein stability value. Variation in Δ*G* among environments was linked to the balanced contribution of hydrophobic amino acids promoting folding*versus* those lowering thermal stability (Table S3).

Although for key enzymes it was found that adaptive differences among species at different depths showed structural changes as well as structural stability (reviews in: \[[@B4]--[@B6]\]), several studies remarked on the importance of regulatory regions of the genome acting on gene regulation (see \[[@B46]\]). Promoter regions often work together with other*cys*-regulatory elements (e.g., transcription factors binding sites, enhancers, silencers, and insulators) to regulate the transcription and expression of mRNA in a specific tissue at different times and places throughout the genome. In addition, in regards to adaptation to depth, it is remarked the importance of osmolyte concentrations of methylamines (TMAO is the most relevant). These are protein stabilizers that counteract inhibition of proteins by hydrostatic pressure \[[@B47]--[@B49]\].

Finally, we should also take into account the stress imposed on this fishes being captured at a depth range of 1100--1250 m and brought to surface in at least 3--5 hours. By the time specimens of black scabbardfish reach the surface they are already dead or dying; therefore the expression of some of the transcripts might be different if compared to a transcriptome of a fish sampled at a thousand meters depth. Notwithstanding, several experiments using model species have been targeting specific genes linked to stress factors \[[@B50]\] and this might prove useful in interpreting transcriptomes of animals undergoing stressful conditions before being sampled.

On the bases of this preliminary but encouraging results, we are planning to explore further deep water adaptivity employing different NGS sequence technology and experimental design, for instance using individual RADseq approach on a larger number of specimens caught at deep or shallow waters or by comparing sister species that live at different depth levels.

3.4. EST-SSR Resources for Population Genetics {#sec3.4}
----------------------------------------------

In this study, a total of 153 EST-SSRs were detected in 142 contigs (3.87%) with a frequency of one EST-SSR per 27,69 kb sequence ([Table 8](#tab8){ref-type="table"}). A further selection was made taking into account the type of repeats, size of the fragment, and quality of the primer sets, reducing the EST-SSR to 98. Some of those microsatellite markers were found in one tissue but not in others (16.3%), while others resolved to be found in more than one tissue (83.7%). Among the identified EST-SSRs, tri-nucleotide repeats represented the largest portion (40.8%), followed by di-nucleotide (28.6%), and tetra-nucleotide (25.7%), and only 2 EST-SSRs penta-nucleotide were identified (Table S4). Primer pairs have sizes ranging between 19 and 25 bp and melting temperatures from 57.1°C to 61.2°C, while expected PCR products are 100 to 280 bp in length. None of these EST-derived primer pairs have been tested for neutrality; therefore not all may be suitable for basic population genetic studies.

4. Conclusions {#sec4}
==============

The transcriptome analysis of*A. carbo*, revealed a comprehensive set of genes expressed in six tissues, producing over 8,000 genes, 55.6% of which annotated by similarity to known proteins or nucleotide sequences from freely accessible database. The transcriptome of black scabbardfish sets the stage for expanding the genetic resources available while providing sets of genes that are likely tied with physiological adaptations to depth, particularly to low temperature and high pressure factors. This study represents the first transcriptome analysis for a deep-sea fish providing insights based on comparison analyses of homologous depth-related functional genes from shallow, deep-water, and abyssal fish highlighting similarities for*A. carbo* isozyme patterns and stability to other bathypelagic fishes. Direct sequence comparison suggested that the signals embedded in these functional genes reflect evolutionary divergence among taxa rather than any kind of enzyme adaptations. Organisms are adapted to deep-sea environment and their physiological tolerance may vary among taxa as well as their enzymatic activities under extreme conditions. Osmolyte concentrations as protein stabilizers that counteract inhibition of proteins by hydrostatic pressure \[[@B6], [@B47]--[@B49]\] play a key role in deep-sea adaptation. Furthermore, contribution to adaptation is also provided by promoter regions that together with*cys*- and*trans*-regulatory elements work concertedly, at the mRNA transcription level, to drive the expression of a gene in a specific tissue or at a specific time \[[@B46]\].

The strong correlation detected between the values of standardized contigs frequency with the expression level of the gene by qPCR supports the use of our data to explore gene expression patterns in*A. carbo*.

Finally, considering the importance of this species for fishery management, a further exploration of this database has provided the characterization of new EST-SSR markers as additional resource for basic population genetic studies as well as tagging and mapping of genes.

Supplementary Material {#supplementary-material-sec}
======================

###### 

As supplementary material we included: amplification plots of all genes selected for this study in different tissue from *Aphanopus carbo* (Figure S1); a summary table of contig frequency in all tissues for those genes associated to depth not listed in Table 3 (Table S1); details on outputs related to the selection models for four teleost genes associated to depth (Table S2); values of Gibbs free energy calculations in terms of kinetic and thermodynamic quantities for all depth-related genes (Table S3); and the list of EST-SSR markers (Table S4).

All molecular work was supported by the research Projects DEECON (European Science Foundation, under the EUROCORES programme, proposal no. 06-EuroDEEP-FP-008) and ReDEco (MarinERA programme funded by the EU FP6 ERA-NET Scheme, MARIN-ERA/MAR/0003/2008). Specimens of*A. carbo*were provided by the CONDOR Project (EEA Grants Financial Mechanism, Iceland, Liechtenstein and Norway, proposal no. PT0040). The authors are thankful to captain, crew, and technicians onboard of the R/V Arquipelago for their excellent contribution while at sea. Sergio Stefanni is a research fellow supported by the Marie Curie grant cofunded by the EU under the FP7-People-2012-COFUND; Cofunding of Regional, National and International Programmes, GA no. 600407 and the Bandiera Project RITMARE. IMAR/DOP is funded through the pluri-annual and programmatic funding scheme as research unit number 531 and associate laboratory number 9. The authors are also grateful to Conceição Egas for the support provided at Biocant. Last, but not least, the authors thank three anonymous reviewers for their comments and suggestions that greatly improved the paper.

Conflict of Interests
=====================

The authors declare that there is no conflict of interests regarding the publication of this paper.

Authors\' Contributions
=======================

Sergio Stefanni conceived the study and design, collected samples, participated in the data analysis, and drafted the paper. Raul Bettencourt contributed to gene validation experiments and discussion. Miguel Pinheiro and Gianluca De Moro developed the pipe-line analysis for all functional annotations and the SQL database. Lucia Bongiorni and Alberto Pallavicini participated in the structure, discussion, and paper drafting.

![Functional categorization of unigenes with gene ontology (GO) term for the*Aphanopus carbo* EST collection. These unigenes results were functionally classified from the six tissues pooled together as percentages under three main functional categories with respective GO Slim terms. Data refer to assemblage derived by implementation of Newbler, Roche 454 sequence analysis software.](IJG2014-267482.001){#fig1}

![Gene expression normalized to a relative value of 1.0 of all genes selected for this study in different tissues from*Aphanopus carbo*. S = spleen, B = brain, H = heart, G = gonad, L = liver, M = muscle,  \*= expression in the brain was below detection and  \*\*= expression in the gonads was below detection.](IJG2014-267482.002){#fig2}

![Comparative plots of relative expression as contig frequencies*versus* mean qPCR values for all genes selected for this study. Codes for each of the gene are listed in [Table 1](#tab1){ref-type="table"}, S = spleen, B = brain, H = heart, G = gonad, L = liver, and M = muscle; the correlation coefficient (Pearson\'s *r*) resulted to be highly significant (*P* \< 0.0001) for most of the genes except for \**P* \< 0.05 and °not significant.](IJG2014-267482.003){#fig3}

![At the top: molecular phylogenetic tree of COI gene estimated by the HKY nucleotide substitution model constructed by neighbor joining and rooted including the sequence of*Latimeria chalumnae* (acc. nr.: NC_001804); in the middle: NJ tree of LDH and MDH genes; and at the bottom: NJ tree of globin (Hb) gene. Numbers in proximity of the nodes denote the bootstrap value (above 50%) out of 1000 replicates. The scale indicates the evolutionary distance of the base substitution per site. Taxon coding includes the first two letters for the genus followed by three letters for the species name (see [Table 6](#tab6){ref-type="table"}). A black square helps to locate*Aphanopus carbo*in the trees.](IJG2014-267482.004){#fig4}

![Comparative plots of normalised difference in protein stability (Δ*G*) for different functional genes among several representative species. Taxon coding includes the first two letters for the genus followed by three letters for the species name (see [Table 6](#tab6){ref-type="table"}).](IJG2014-267482.005){#fig5}

###### 

List of targeted genes using qPCR, primer sets specifically designed for this study, size for each of the product, and NCBI accession number for all the EST sequences.

  Gene                                                       Primer name               Primer Sequences (5′-3′)   Size (bp)   NCBI Accession \#
  ---------------------------------------------------------- ------------------------- -------------------------- ----------- --------------------
  Elongation factor 1-beta                                   *EF-1B L*                 GCTTGGACATGTCGGTCTCGTC     229 bp      All_gs454_000396
  *EF-1B H*                                                  GTGGCTGACACCACATCTGGC                                            
  Ras-related GTP-binding protein A                          *Rab-1A L*                AGTAGCCGTTCCACCTTGTCGG     247 bp      All_gs454_000598
  *Rab-1A H*                                                 TGCCAAGAAACCGTACGTGGGA                                           
  Basic Transcription Factor 3-like 4                        *BTF3 L*                  CCCAAAGTTCAGGCCTCCCTGT     273 bp      All_gs454_000873
  *BTF3 H*                                                   TCATGTGCGTCAGTTCGCTTCG                                           
  Cu/Zn Superoxide Dismutase                                 *SOD-1 L*                 AAACGTGACTGCAGGAGGGGAT     240 bp      All_gs454_000925
  *SOD-1 H*                                                  CAGTGCTCCTGCTCCATGTTCG                                           
  2-Cys Peroxiredoxin                                        *PRDX1 L*                 CCGATAACCTCGCAGCCGATAC     243 bp      All_gs454_000558
  *PRDX1 H*                                                  ACAGTCATTTGCCACCAGCATCA                                          
  Heat Shock Protein 90                                      *HSP90 L*                 TGACGATGTCCCCACAGATGAGG    221 bp      All_gs454_000008
  *HSP90 H*                                                  GCAACACTGGTCCACCACACAAC                                          
  Ferritin, heavy subunit                                    *Ferr L*                  CCTGCAGCTTGAGAAGAGCGTC     203 bp      All_gs454_000681
  *Ferr H*                                                   CAAACAGGTACTCGGCCATGCC                                           
  *α* ~2~ Globin                                             *Hb-A L*                  AAATTGTTGGCCATGCGGAGGA     208 bp      All_gs454_001919
  *Hb-A H*                                                   CTGAGGTTCAGCAGACCTGCCT                                           
  *β* ~2~ Globin                                             *Hb-B L*                  TCGTCTACCCCTGGTGTCAGAG     245 bp      All_gs454_001018
  *Hb-B H*                                                   AACCACAATGGTCAGGCAGTCC                                           
  Ependymin-1 precursor                                      *EPD-1 L*                 CAGGTGTGAGGCAGTGCAGT       230 bp      All_gs454_000469
  *EPD-1 H*                                                  ACCCCGATCTCCTCCTGGTG                                             
  Fatty acid-binding protein, brain                          *BLBP L*                  CAACACTTCTTGGCCGGTTTGG     239 bp      All_gs454_001220
  *BLBP H*                                                   GAGAGGAGTTCGACGAAGCCAC                                           
  CD63-like protein Sm-TSP-2                                 *TSPAN-8 L*               TCGCTGGCTGCTCTGAGAAAGA     200 bp      All_gs454_000381
  *TSPAN-8 H*                                                GGTCACGCCGAGCTGTATTCTG                                           
  Tropomyosin 4 isoform 1                                    *TRPM-1 L*                GTGGAGGAGGAGTTGGACCGAG     221 bp      All_gs454_000222
  *TRPM1 H*                                                  TTGCGAGCCACCTCCTCGTATT                                           
  C-Myc-binding protein                                      *MYCBP L*                 CGCCAGTTTACCTGCGTTCCAA     182 bp      All_gs454_001640
  *MYCBP H*                                                  GGCCGTCAACAACACCACCTTT                                           
  Cathepsin S                                                *CTSS L*                  AACAGCCTACCCCTACACAGCC     200 bp      All_gs454_000156
  *CTSS H*                                                   TGTACACACCGTGGCGGTAGAA                                           
  Transferrin                                                *STF-1 L*                 AGCTGCACCAGCTTCACAGTTG     215 bp      All_gs454_000004
  *STF-1 H*                                                  AAGGATGGCACCAGACAACCCA                                           
  Warm Temperature Acclimation related-like 65 kDa protein   *HPX L*                   TGATACCGGGTGGAACCTGGTG     207 bp      All_gs454_000060
  *HPX H*                                                    GCTGCTGTGGAGTGTCCCAAAG                                           
  Betaine Homocysteine S-methyltransferase                   *BHTM L*                  GGGGGTTCGCTGTTACCAAGTG     194 bp      All_gs454_000088
  *BHTM H*                                                   TGTGAGACAGCAGCCTCAGGAG                                           
  FUCL1 Fucolectin-1                                         *FUCL1 L*                 CGCAAACCCTTTGGCTGGTGTA     196 bp      All_gs454_000758
  *FUCL1 H*                                                  GGCTTTTCCTTGGACTGCCAGG                                           
  Aldolase B                                                 *ALDB L*                  GCCATTGGTCTTGGCCCTGATC     220 bp      All_gs454_000115
  *ALDB H*                                                   CGCTGTGCCTGGTATCTGCTTC                                           
  Type-4 ice-structuring protein LS-12 precursor             *ISP LS12 L*              AAGACCTGACAAACCAGGCCCA     198 bp      All_gs454_001277
  *ISP LS12 H*                                               GGAGGATGGCCTCCATCTGCTT                                           
  Alcohol Dehydrogenase 8a                                   *ADH L*                   GGCAAGAAGGTGCTGCAGTTCA     228 bp      All_gs454_000105-6
  *ADH H*                                                    CATGACTGCAGCCAAACCCACA                                           
  Glyceraldehyde-3-phosphate Dehydrogenase                   *GAPDH L*                 GTCAACCACTGACACGTTGGGG     229 bp      All_gs454_000148
  *GAPDH H*                                                  CGGCATCATTGAGGGCCTGATG                                           
  Lactate Dehydrogenase-A                                    *LDH-A L*                 TCTTAACCTGGTGCAGCGCAAC     219 bp      All_gs454_000149
  *LDH-A H*                                                  TGGAGCTTCTCGCCCATGATGT                                           
  Phosphoglycerate Mutase 2-1 (muscle)                       *PglyM L*                 ACACCTCTGTGCTGAAACGTGC     212 bp      All_gs454_000309
  *PglyM H*                                                  CATGGGTGGAGGTGGGATGTCA                                           
  Heat Shock Protein 70                                      *HSP70 L*                 CGGTGTTGTGTGCTGGGTGAAA     207 bp      All_gs454_000005
  *HSP70 H*                                                  CCACATAGCTGGGTGTGGTCCT                                           
  Fructose-bisphosphate Aldolase A                           *FBPA L*                  GGAACCAACGGCGAGACAACAA     208 bp      All_gs454_002732
  *FBPA H*                                                   CAATGGGGACGATGCCATGCAT                                           
  Phosphoglucose Isomerase-2                                 *PGI L*                   CCACACTGGGCCAATTGTCTGG     217 bp      All_gs454_000011
  *PGI H*                                                    GGCCTCCTCTGTGGTCTTACCC                                           

###### 

Global statistics for each of the nonnormalized libraries using Newbler software.

  Tissue                        Spleen      Brain       Heart        Gonad        Liver        Muscle       Total
  ----------------------------- ----------- ----------- ------------ ------------ ------------ ------------ -------------
  Total EST                     15,034      33,337      73,263       157,275      134,523      92,788       544,491
  Total bases                   3,426,510   8,219,500   17,647,600   37,123,700   31,792,600   23,342,900   129,412,000
  Contigs                       567         651         1,260        3,875        1,274        626          8,319
  Average contig length         470         619         567          465          612          689          555
  Contigs *e*-value \< 10^−6^   220         420         622          951          617          409          2,440
  ESTscan                       584         345         977          2,838        634          269          2,715
  No similarity                 36          70          109          406          211          74           1,128
  GO annotation                 202         338         473          623          509          307          1,728
  InterPro annotation           223         417         610          908          649          395          2,395

###### 

Similarity search for unigenes comparing transcriptomes of other teleosts against *Aphanopus carbo* using two levels of stringency.

  ---------------------------------------------------------------------------------------------
  Species                    Sequences available   *A. carbo*  \         *A. carbo*  \
                                                   *e*-value \< 10^−3^   *e*-value \< 10^−10^
  -------------------------- --------------------- --------------------- ----------------------
  *Danio rerio*              42,787                2,457                 2,199

  *Gasterosteus aculeatus*   27,576                2,445                 2,223

  *Oreochromis niloticus*    26,763                2,436                 2,202

  *Oryzias latipes*          24,674                2,412                 2,173

  *Takifugu rubripes*        47,841                2,268                 2,062

  *Tetraodon nigroviridis*   23,118                2,300                 2,073
  ---------------------------------------------------------------------------------------------

###### 

List of Pfam conserved domains that were tissue specifically expressed in *Aphanopus carbo* transcriptome. Frequencies and *P*-values of tissue specificity analysis are indicated.

  Gonads    Heart            Liver                                                                                             
  --------- ---------------- ------- -------------- --------- --------------- --- ------------- --------- --------------- ---- --------------
  PF00100   Zona pellucida   17      3.15 10^−13^   PF00011   HSP20           3   4.03 10^−4^   PF00084   Sushi           14   1.44 10^−10^
  PF00125   Histone          8       6.06 10^−5^    PF13405   EF hand 4       5   6.77 10^−4^   PF00089   Trypsin         20   7.59 10^−9^
  PF01400   Astacin          3       3.22 10^−3^    PF00412   LIM             3   1.52 10^−3^   PF00079   Serpin          12   5.63 10^−6^
  PF00069   Pkinase          3       0.01           PF05556   Calsarcin       3   3.60 10^−3^   PF07678   A2M comp        4    1.35 10^−4^
  PF00653   BIR              2       0.02           PF01576   Myosin tail 1   3   3.60 10^−3^   PF01042   Ribonuc L-PSP   3    1.26 10^−3^
  PF13424   TPR 12           2       0.02           PF00056   Ldh 1 N         3   3.60 10^−3^   PF00045   Hemopexin       3    1.26 10^−3^
  PF13695   zf-3CxxC         2       0.02           PF05300   DUF737          2   5.50 10^−3^   PF00059   Lectin C        6    1.69 10^−3^
  PF09360   zf-CDGSH         2       0.02           PF00992   Troponin        4   6.40 10^−3^   PF00701   DHDPS           4    4.64 10^−3^
  PF01712   dNK              2       0.02           PF00595   PDZ             3   6.82 10^−3^   PF00386   C1q             8    6.63 10^−3^
  PF00538   Linker histone   2       0.02           PF00022   Actin           3   0.01          PF00021   UPAR LY6        5    0.01
  PF10178   DUF2372          2       0.02           PF02874   ATP-synt ab N   2   0.02          PF00754   F5 F8 type C    9    0.01
  PF04856   Securin          2       0.02           PF13895   Ig 2            3   0.02          PF08702   Fib alpha       2    0.01
  PF00250   Fork head        2       0.02           PF00191   Annexin         2   0.03          PF03982   DAGAT           2    0.01
  PF01498   HTH Tnp Tc3 2    3       0.03           PF00212   ANP             2   0.03          PF01048   PNP UDP 1       2    0.01
  PF00268   Ribonuc red sm   2       0.06           PF05347   Complex1 LYR    2   0.07          PF01014   Uricase         2    0.01

  Muscle    Brain             Spleen                                                                                            
  --------- ----------------- -------- ------------- --------- ---------------- --- ------------- --------- --------------- --- -------------
  PF01410   Troponin          7        1.99 10^−5^   PF01669   Myelin MBP       4   4.03 10^−5^   PF00042   Globin          5   1.33 10^−6^
  PF02807   COLFI             3        9.67 10^−4^   PF01453   B lectin         2   6.44 10^−3^   PF00078   RVT 1           5   1.33 10^−6^
  PF00041   ATP-gua PtransN   3        3.59 10^−3^   PF00612   IQ               2   6.44 10^−3^   PF00993   MHC II alpha    4   5.02 10^−6^
  PF02453   fn3               3        3.59 10^−3^   PF11414   Suppressor APC   2   6.44 10^−3^   PF07686   V set           5   2.49 10^−6^
  PF13895   Reticulon         3        3.59 10^−3^   PF05768   DUF836           2   6.44 10^−3^   PF07654   C1 set          5   4.71 10^−4^
  PF00365   Ig_2              4        7.97 10^−3^   PF05196   PTN MK N         2   6.44 10^−3^   PF00089   Trypsin         5   2.01 10^−3^
  PF01216   PFK               2        9.84 10^−3^   PF04300   FBA              2   6.44 10^−3^   PF01391   Collagen        2   2.30 10^−3^
  PF01267   Calsequestrin     2        9.84 10^−3^   PF00287   Na K-ATPase      2   6.44 10^−3^   PF00240   ubiquitin       3   3.28 10^−3^
  PF00261   F-actin cap A     2        9.84 10^−3^   PF11032   ApoM             3   8.53 10^−3^   PF09307   MHC2-interact   2   6.67 10^−3^
  PF00856   Tropomyosin       2        9.84 10^−3^   PF00061   Lipocalin        4   0.01          PF00643   Zf-B box        2   0.01
  PF01661   SET               2        9.84 10^−3^   PF00007   Cys knot         2   0.02          PF13445   Zf-RING LisH    2   0.01
  PF01667   Macro             2        0.03          PF01275   Myelin PLP       2   0.02          PF01498   HTH Tnp Tc3 2   2   0.02
  PF00036   Ribosomal S27e    2        0.05          PF00300   His Phos 1       2   0.02          PF02301   HORMA           1   0.05
  PF01576   efhand            2        0.05          PF00230   MIP              2   0.02          PF14259   RRM 6           1   0.05
  PF01410   Myosin_tail_1     2        0.08          PF00091   Tubulin          3   0.02          PF09004   DUF1891         1   0.05

###### 

Summary table of contig frequency for each of the tissues for those genes used to test the correlation with qPCR assays. S: spleen, B: brain, H: heart, G: gonad, L: liver and M: muscle. Names of genes based on their coding used in this table are found in table.

  Gene       Contig code       Spleen   Brain   Gonads   Heart   Liver   Muscle
  ---------- ----------------- -------- ------- -------- ------- ------- --------
  EF-1B      isotig00991       9        24      33       19      38      54
  Rab-1A     isotig02689       4        0       3        3       0       4
  BTF3       isotig02213       3        5       7        10      8       2
  SOD-1      isotig02222       4        16      32       16      25      14
  PRDX1      isotig01838       4        50      59       17      23      27
  HSP90      isotig01406       3        64      94       86      86      13
  Ferr       isotig01665       33       121     10       210     264     47
  Hb-A       isotig06973       406      17      2        24      108     1
  Hb-B       isotig00163       2,303    81      16       313     707     22
  EPD-1      isotig01567       0        1,190   0        1       0       0
  BLBP       isotig02632       0        171     0        0       0       0
  TSPAN-8    isotig01397       17       7       3        407     24      9
  TRPM-1     isotig01595       1        4       0        637     2       2
  MYCBP      isotig01988       0        2       135      1       1       1
  CTSS       isotig01659       0        0       135      0       0       0
  STF-1      isotig00767       0        1       30       0       3,218   0
  HPX        isotig01479       0        0       0        0       1,234   0
  BHTM       isotig01489       1        0       0        0       1,024   0
  FUCL1      isotig00473       1        0       0        0       22      0
  ALDB       isotig01524       0        1       30       0       369     0
  ISP LS12   isotig03004       0        0       0        0       211     1
  ADH        isotig01491-546   1        2       16       7       292     6
  GAPDH      isotig00609       0        8       51       539     134     1,281
  LDH-A      isotig01401       2        7       2        5       0       789
  PglyM      isotig01642       0        0       0        36      0       241
  HSP70      isotig01398       0        0       0        0       2       142
  FBPA       isotig00492       2        2       0        233     0       4,471
  PGI        isotig01410       0        0       0        26      0       151

###### 

Species list with informations on their type of environment (M = marine, BR = brackish and FW = freshwater), climate, depth range (in meters) and the type of genes used in the study with relative Genbank accession numbers.

  Order                Family            Species                        Common name                Environment   Climate       Depth range         Gene                       NCBI Acc. Nr
  -------------------- ----------------- ------------------------------ -------------------------- ------------- ------------- ------------------- -------------------------- ----------------------------------------------------
  Perciformes          Trichiuridae      *Aphanopus carbo*              Black scabbardfish         M             Deep-water    200--1700           COI                        EU854076
  Beloniformes         Adrianichthydae   *Oryzias latipes*              Japanese rice fish         FW + BR       Subtropical   shallow             ACTA1, MDHc, MyHC, COI     NM_001104806, NM_001163134, XM_004071618, AB498066
  Scorpaeniformes      Hexagrammidae     *Pleurogrammus azonus*         Okhotsk atka mackerel      M             Temperate     0--240              ACTA1                      AB073381
  Perciformes          Scombridae        *Scomber scombrus*             Atlantic mackerl           M             Temperate     0--200 (0--1000)    ACTA1, COI                 EF607093, KC015895
  Perciformes          Percihcthyidae    *Siniperca chuatsi*            Mandarin fish              FW            Temperate     10                  ACTA1, MyHC, COI           AY395872, AY454304, NC_015822
  Perciformes          Sparidae          *Sparus aurata*                Gilthead seabream          M             Temperate     1--30 (1--150)      ACTA1                      AF190473
  Perciformes          Sphyraenidae      *Sphyraena idiastes*           Pelican barracuda          M             Tropical      3--24               ACTA1, LDH-A, MDHc, mMDH   AF503593, SIU80001, AF390559, AF390561
  Tetraodontiformes    Tetraodontidae    *Takifugu rubripes*            Japanese pufferfish        M + FW + BR   Temperate     0--200 (0--1000)    Hb-A, mMDH, COI            XM_003964767, XM_003965959, HM102315
  Scorpaeniformes      Anoplopomatidae   *Anoplopoma fimbria*           Sablefish                  M             Deep-water    0--2740             Hb-B, COI                  BT082849, JQ353978
  Perciformes          Serranidae        *Epinephelus coioides*         Orange-spotted grouper     M + BR        Subtropical   1--100              Hb-B                       GU982530
  Gasterosteiformes    Gasterosteidae    *Gasterosteus aculeatus*       Three-spined stickleback   M + FW + BR   Temperate     0--100              Hb-B                       NM_001267638
  Perciformes          Nototheniidae     *Notothenia coriiceps*         Black rockcod              M             Polar         0--550              LDH-A, MyHC, COI           AF079822, AJ243767, EU326390
  Perciformes          Nototheniidae     *Notothenia angustata*         Maori chief                M             Temperate     0--100              Hb-A, Hb-B                 P62363, P29628
  Gadiformes           Macrouridae       *Coryphaenoides armatus*       Abyssal grenadier          M             Deep-water    282--5180           LDH-B, MyHC, COI           AJ609232, AB330140, FJ164497
  Gadiformes           Gadidae           *Gadus morhua*                 Atlantic cod               M + BR        Temperate     0--600              LDH-B, COI                 AJ609233, KC015385
  Gadiformes           Gadidae           *Arctogadus glacialis*         Arctic cod                 M             Deep-water    0--1000             Hb-A, COI                  Q1AGS4, KC015200
  Perciformes          Latidae           *Lates calcarifer*             Barramundi                 M + FW + BR   Tropical      10--40              LDH-B, COI                 FJ439507, JQ431879
  Gadiformes           Gadidae           *Merlangius merlangus*         Whiting                    M             Temperate     30--100 (10--200)   LDH-B, COI                 AJ609234, JQ623954
  Cyprinodontiformes   Poeciliidae       *Poecilia reticulata*          Guppy                      FW + BR       Tropical      Shallow             LDH-B, COI                 EF408825, JX968696
  Gadiformes           Macrouridae       *Trachyrincus murrayi*         Roughnose grenadier        M             Deep-water    0--1630             LDH-B, COI                 AJ609235, AP008990
  Perciformes          Channichthyidae   *Chionodraco rastrospinosus*   Ocellated icefish          M             Polar         0--1000             LDH-A, COI                 AF079829, EU326337
  Perciformes          Pomacentridae     *Chromis caudalis*             Blue-axil chromis          M             Tropical      15--55              LDH-A                      AY289558
  Perciformes          Gobiidae          *Rhinogobiops nicholsii*       Blackeye goby              M             Subtropical   ?-106               LDH-A                      AF079534
  Cypriniformes        Cyprinidae        *Cyprinus carpio*              Common carp                FW + BR       Subtropical   shallow             LDH-A, MyHC, COI           AF076528, D89992, HQ960709
  Cyprinodontiformes   Fundulidae        *Fundulus heteroclitus*        Mummichog                  M + FW + BR   Temperate     shallow             LDH-A, LDH-B, COI          L43525, L23792, EU524629
  Osmeriformes         Osmeridae         *Osmerus mordax*               Rainbow smelt              M + FW + BR   Temperate     0--425              MDHc, mMDH                 BT075651, BT075600
  Salmoniformes        Salmonidae        *Salmo salar*                  Atlantic salmon            M + FW + BR   Temperate     10--23 (0--210)     MDHc, mMDH                 BT060183, BT048216
  Gadiformes           Macrouridae       *Coryphaenoides acrolepis*     Pacific grenadier          M             Deep-water    900--1300           MyHC, COI                  AB330141, JQ354060
  Gadiformes           Macrouridae       *Coryphaenoides yaquinae*      n.a.                       M             Deep-water    3400--5800          MyHC, COI                  AB330139, GU440291
  Perciformes          Cirrhitidae       *Paracirrhites forsteri*       Blackside hawkfish         M             Tropical      5--35               MyHC, COI                  AJ243770, HQ561521
  Perciformes          Carangidae        *Seriola dumerili*             Greater amberjack          M             Subtropical   18--72 (1--360)     MyHC, COI                  AB032020, KC015917
  Gadiformes           Gadidae           *Boreogadus saida*             Polar cod                  M + BR        Polar         0--400              Hb-A, Hb-B, COI            DQ125471, Q1AGS6, KC015250

###### 

Comparison of protein sequences of depth related genes between *Aphanopus carbo* database (All_gs454_xxx) and other fishes. Diff.: amino acid differences; Id. %: percentage of identity; Gaps: number of gaps introduced in the complete alignment and Acc. nr.: NCBI Accession number. ^1^Partial sequence; ^2^only AA sequence available.

###### 

\(a\) LDH-A (332 AA)

  All_gs454_00149*vs.*           Diff.   Id. %   Gaps   Acc. nr.
  ------------------------------ ------- ------- ------ ----------
  *Sphyraena idiastes*           16      95.18   0      U80001
  *Rhinogobiops nicholsii*       17      94.88   0      AF079534
  *Chromis caudalis*             21      93.67   0      AY289558
  *Chionodraco rastrospinosus*   26      92.17   1      AF079829
  *Notothenia coriiceps*         27      91.87   1      AF079822
  *Fundulus heteroclitus*        27      91.87   0      L43525
  *Cyprinus carpio*              44      86.75   0      AF076528

###### 

\(b\) MDHc-A (333 AA)

  All_gs454_00146*vs.*   Diff.   Id. %   Gaps   Acc. nr.
  ---------------------- ------- ------- ------ ------------
  *Oryzias latipes*      15      95.50   0      NM_1163134
  *Sphyraena idiastes*   20      93.99   0      AF390559
  *Osmerus mordax*       30      90.99   0      BT075651
  *Salmo salar*          35      89.49   0      BT060183

###### 

\(c\) Actin-1 (375 AA)

  All_gs454_00104*vs.*            Diff.   Id. %   Gaps   Acc. nr.
  ------------------------------- ------- ------- ------ -----------------
  *Scomber scombrus 1*            0       100     0      EF607093
  *Iniperca chuatsi 1*            0       100     0      AY395872
  *Oryzias latipes 1*             1       99.73   0      NM_1104806
  *Pleurogrammus azonus 1*        1       99.73   0      AB073381
  *Coryphaenoides acrolepis 1*    1       99.73   0      AB021649
  *Coryphaenoides cinereus 1*     1       99.73   0      AB021651
  *Coryphaenoides armatus 2a*     2       99.47   0      AB086240
  *Coryphaenoides yaquinae 2a*    2       99.47   0      AB086242
  *Coryphaenoides acrolepis 2a*   2       99.47   0      AB021650
  *Coryphaenoides cinereus 2a*    2       99.47   0      AB021652
  *Cyprinus carpio 1*             2       99.47   0      AY395870
  *Sphyraena idiastes 1*          3       99.20   0      AF503593
  *Coryphaenoides armatus 2b*     5       98.67   0      AB086241
  *Coryphaenoides yaquinae 2b*    5       98.67   0      AB086243
  *Aphanopus carbo 2a*            6       98.40   0      All_gs454_00102
  *Sparus aurata 1*               9       97.60   0      AF190473

###### 

\(d\) LDH-B (334 AA)

  All_gs454_00143*vs.*       Diff.   Id. %   Gaps   Acc. nr.
  -------------------------- ------- ------- ------ ----------
  *Lates calcarifer*         14      95.81   0      FJ439507
  *Fundulus heteroclitus*    21      93.71   0      L23792
  *Trachyrincus murrayi*     47      85.93   0      AJ609235
  *Coryphaenoides armatus*   50      85.03   0      AJ609232
  *Merlangius merlangus*     52      84.43   1      AJ609234
  *Gadus morhua*             55      83.53   1      AJ609233

###### 

\(e\) MyHC^1^ (305 AA)

  All_gs454_00001*vs.*         Diff.   Id. %   Gaps   Acc. nr.
  ---------------------------- ------- ------- ------ ------------
  *Paracirrhites forsteri*     53      94.95   0      AJ243770
  *Seriola dumerilii*          54      94.86   0      AB032020
  *Siniperca chuatsi*          59      94.38   0      AY454304
  *Oryzias latipes*            62      94.10   2      XM_4071618
  *Cyprinus carpio*            69      93.43   2      D89992
  *Coryphaenoides acrolepis*   86      91.82   1      AB330141
  *Notothenia coriiceps*       86      91.82   1      AJ243767
  *Coryphaenoides yaquinae*    89      91.53   1      AB330139
  *Coryphaenoides armatus*     92      91.25   1      AB330140

###### 

\(f\) Hb-A^1^ (143 AA)

  All_gs454_01074*vs.*     Diff.   Id. %   Gaps   Acc. nr.
  ------------------------ ------- ------- ------ ------------
  *Notothenia angustata*   39      72.73   1      P62363^2^
  *Boreogadus saida*       43      69.93   0      DQ125471
  *Arctogadus glacialis*   44      69.93   0      DQ125475
  *Takifugu rubripes*      46      67.83   0      XM_3964767
  *Gadus morhua*           53      62.94   0      O42425^2^

###### 

\(g\) Hb-B^1^ (146 AA)

  All_gs454_01018*vs.*       Diff.   Id. %   Gaps   Acc. nr.
  -------------------------- ------- ------- ------ ------------
  *Epinephelus coioides*     27      81.63   0      GU982530
  *Anoplopoma fimbria*       31      78.91   0      BT082849
  *Gasterosteus aculeatus*   31      78.91   0      NM_1267638
  *Boreogadus saida*         40      72.79   0      Q1AGS6^2^
  *Notothenia angustata*     42      71.43   0      P29628^2^

###### 

Statistics of EST-SSRs identified in *Aphanopus carbo* transcriptome.

  Searching item                                   Numbers
  ------------------------------------------------ -----------
  Total number of sequences examined               7,920
  Total size of examined sequences (bp)            4,235,839
  Total number of identified SSRs                  153
  Number of SSR containing sequences               142
  Number of sequences containing more than 1 SSR   9
  Number of SSRs present in compound formation     8
  Di-nucleotide                                    63
  Tri-nucleotide                                   52
  Tetra-nucleotide                                 35
  Penta-nucleotide                                 3

[^1]: Academic Editor: Elena Pasyukova
