Procedures for fish manipulation and tissue collection were approved by the Ethics and Animal Welfare Committee of Institute of Aquaculture Torre de la Sal and carried out according to the National (Royal Decree RD53/2013) and the current EU legislation (2010/63/EU) on the handling of experimental fish.

Fish were reared from early life stages under natural conditions of photoperiod and temperature at the experimental facilities of IATS (40°5 N; 0°10 E). Blood of one single male was obtained from caudal vessels using heparinized syringes, and DNA from total blood cells was extracted with a commercial kit (RealPure Spin Blood Kit, Durviz, Valencia, Spain). Quality and quantity of genomic DNA was assessed by means of PicoGreen quantification and gel electrophoresis. An aliquot of 5 μg DNA was mechanically sheared with a bath sonicator (Diagenode BioRuptor, Diagenode, Liège, Belgium) and low molecular weight fragments were used for the preparation of DNA libraries.

Total RNA (70–100 μg) from white skeletal muscle (six individual juvenile fish with a body weight of 78.6 ± 5.2 g) and samples of anterior and posterior intestine sections was extracted with the MagMAX^TM-96 Total RNA Isolation Kit (Applied Biosystems, Foster City, CA, United States). The RNA concentration and purity was determined using a Nanodrop 2000c (Thermo Scientific, Wilmington, DE, United States). Quality and integrity of the isolated RNA were checked on an Agilent Bioanalyzer 2100 total RNA Nano series II chip (Agilent, Amstelveen, Netherlands), yielding RNA integrity numbers (RINs) between 8 and 10.

Genomic DNA material was used for the preparation of two standard TrueSeq Illumina libraries (Illumina, Inc.) with an average size of 360 and 747 bp, respectively. Illumina NextSeq500 system under a 2 × 150 paired-end (PE) format was used as sequencing platform to generate approximately 600 million reads. Additionally, two different strategies were implemented in order to help in genome scaffolding: (1) Nextera Mate-Pair Preparation Kit (Illumina, Inc.) was used to make two mate pairs (MP) libraries (average insert sizes were 5 and 8 kb) using the Illumina NextSeq500 platform to a depth of 11 Gb (2 × 75 MP format) and (2) genomic DNA was submitted to Macrogen (Seoul, South Korea) for the construction of 12 single molecule real time (SMRT) cell libraries (insert size up to 50 kb) using PacBio RS II (Pacific Biosciences) as sequencing system. Additionally, eight RNA-seq libraries (for more details, see section “Data Availability Statement”) were constructed by means of Illumina TrueSeq RNA-seq preparation protocol (non-directional method). Sequencing of indexed libraries was performed on the Illumina Hiseq v3, resulting in approximately 11–17 million reads per sample (1 × 75 nt single reads) from six skeletal muscle samples and 22–27 million read pairs (2 × 150 nt paired reads) from two pooled intestine samples.

The SMRT cell libraries were pre-processed using the trimming of the CANU assembler (Koren et al., 2017). Illumina PE libraries were checked for quality analysis using FASTQC 0.11.7, available at http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/, and then pre-processed using Cutadapt v1.16 (Martin, 2011) and Prinseq 0.20.4 (Schmieder and Edwards, 2011). Quality analysis and pre-processing of Illumina MP libraries was performed with FastQC and Platanus (Kajitani et al., 2014). These protocols for pre-processing de novo assembly were executed using the DeNovoSeq pipeline provided by the GPRO suite (Futami et al., 2011). Jellyfish (Marçais and Kingsford, 2011) was used to estimate the genome size calculating the count distribution of k-mers in the set of Illumina PE libraries. The estimated coverage was inferred using Bowtie2 v2.3.4.1 (Langmead et al., 2009). Illumina PE and MP libraries were introduced in the 127mer version of the assembler SOAP de Novo2 v2.04-r241 (Luo et al., 2012) for the assembly of gilthead sea bream genome. In order to test different k-mer values, different assemblies were performed and a k-mer length of 63 bp (k63) was considered the best in terms of metrics. To improve the consensus sequence and to close gaps, two rounds of the following combined strategy were conducted: (1) elimination of duplicates with Dedupe of BBTools¹, (2) gap filling using PacBio corrected reads with PBJelly (English et al., 2012), (3) gap filling using PE and MP libraries with Soap de novo Gap Closer, (4) hybrid re-scaffolding using corrected SMRT reads together with Illumina PE and MP reads with Opera 2.0.6 (Gao et al., 2011), and (5) transcriptome-guided re-scaffolding using as reference the gilthead sea bream transcriptome (Calduch-Giner et al., 2013) with L RNA scaffolder (Xue et al., 2013). A step of genome masking was not considered in order to achieve a more reliable genome draft.

Highly conserved non-coding elements (CNEs) present in 3 hermaphrodite genomes (S. aurata, Lates calcarifer, Monopterus albus) were released by Pauletto et al. (2018), and the super-scaffold coordinates related to these CNEs (200–800 bp interval length) were then retrieved. Sequences were aligned against our assembly for increasing the super-scaffolding by means of the BLAST package. A genome browser was built for the navigation and blast-query of the assembled sequences and associated annotations using Javascript-based tool JBrowse (Skinner et al., 2009). The genome browser, available online at www.nutrigroup-iats.org/seabreamdb, provides two modes of navigation for the assembly scaffolds and the entire set of super-scaffolds anchored from CNEs.

Prediction of coding genes was carried out using the software AUGUSTUS 3.3 in a two-step process. An initial round of prediction was conducted, and gene model parameters were trained from a set of 13 fish species (Astyanax mexicanus, Danio rerio, Gadus morhua, Gasterosteus aculeatus, Latimeria chalumnae, Lepisosteus oculatus, Oreochromis niloticus, Oryzias latipes, Petromyzon marinus, Poecilia formosa, Takifugu rubripes, Tetraodon nigroviridis, and Xiphophorus maculatus) available in the Ensembl database release 87 (Cunningham et al., 2015). Then, the merged prediction of gilthead sea bream genes was translated to peptides using OrfPredictor script (Min et al., 2005), and it was used by Scipio 1.4 (Keller et al., 2008) to generate a new training set for a second round of gene predictions. This second round included sequences from the published gilthead sea bream transcriptome (Calduch-Giner et al., 2013) and RNA-seq data from muscle and intestine in addition to those of liver, gills and spleen, retrieved from the SRA archive (see section “Data Availability Statement”) (Piazzon et al., 2019) as AUGUSTUS hints. The script autoAugTrain.pl of AUGUSTUS was used to determine the precise exon/intron gene structures. The Gffread software (Trapnell et al., 2012) rendered the final set of coding sequences (CDSs), using the genome transcript file generated by AUGUSTUS. BLAST package was used for gene annotation, performing BLASTX searches against SWISSPROT, NR and the IATS-CSIC gilthead sea bream transcriptome databases with an E-value cutoff of 10^–5 using the DeNovoSeq pipeline provided by the GPRO suite. Redundancy analysis were performed in order to detect segmental duplications (i.e., predicted genes that occur at more than one site within the genome and typically share > 90% of sequence identity) within the final set of transcripts retrieved from RNA-seq libraries using Dedupe of BBTools¹. Identity thresholds in redundancy analysis were fixed at 90, 95, and 98%.

The mobilome draft was annotated considering the following mobile genetic elements (MGEs): non-coding RNA genes, introns, low complexity repeats, Class I retrotransposons, Class II DNA transposons, and Chimeric/Composite genes. Introns were retrieved from the ab initio predictions. To annotate non-coding RNAs (ncRNAs), a non-redundant database of both small and long ncRNAs was constructed based on the ncRNAs annotations of fish genomes used for de novo gene prediction (Maere et al., 2005). An additional fish tRNA database was created using the tRNAs from D. rerio, G. aculeatus, O. latipes, P. marinus, T. rubripes, and T. nigroviridis from UCSC². Then, a BLAT search (Kent, 2002) served to annotate ncRNAs in the gilthead sea bream genome. Duplicated BLAT outputs were removed using Bedtools³. A final step of curation was performed based on the merging of entries that in the same scaffold had: (1) the same parent and were consecutive in 5–10 nucleotides, (2) the same target and initial position), (3) the same biotype and overlapped, and (4) the support of real transcripts from the gilthead sea bream transcriptome. After curation, repeat sequences retained into longer ones were discarded. To annotate the remaining MGEs, RepeatModeler 1.0.11⁴ was used for the de novo repeat family identification. RepeatMasker 4.0.7 and NCBI-BLAST alignments (E-value threshold < 10^–5) (Altschul et al., 1990) were used to identify simple repeats, low complexity repeats and interspersed repeats within the gilthead sea bream genome. Repbase 22.09 (Bao et al., 2015), GyDB (Llorens et al., 2011), and de novo repeat families coming from RepeatModeler were used as libraries. LTR finder (Xu and Wang, 2007) and Einverted of EMBOSS (Rice et al., 2000) were used to characterize long terminal repeats (LTRs) and inverted repeats, respectively.

All the annotations corresponding to coding genes associated to MGEs (chimeric/composite genes) were extracted from the previously presented coding genes annotation and were used as queries in a BLAST search against Repbase 22.09 and GyDB databases. All the results were curated by means of merging overlapping features with the same annotation or separated by less than 100 nucleotides.

Synteny detection was performed across the genome of gilthead sea bream over other 9 fish species (Cynoglossus semilaevis, D. rerio, G. aculeatus, Maylandia zebra, O. mykiss, O. niloticus, O. latipes, Salmo salar, and X. maculatus). The algorithm includes the following steps: (1) selection of single-copy genes per scaffold in the gilthead sea bream assembly, (2) alignment of gilthead sea bream genes against the other species with BLASTX of the NCBI-BLAST package with more than 70% of sequence identity and coverage, and (3) synteny file construction, establishing an E-value < 10^–5 to consider a gilthead sea bream-species gene correspondence (with number of gaps < 25). A syntenic block must contain a minimum of five genes to be included in the results. Circular genome representations were created using Circos (Krzywinski et al., 2009).

The gilthead sea bream phylome was reconstructed using phylomeDB pipeline (Huerta-Cepas et al., 2014). For each protein-coding gene in gilthead sea bream, a Smith-Waterman search was performed against the proteome database of 19 selected species (S. aurata, Latimeria chalumnae, L. oculatus, D. rerio, A. mexicanus, P. formosa, G. morhua, O. mykiss, Scophthalmus maximus, O. latipes, O. niloticus, T. rubripes, G. aculeatus, T. nigroviridis, Petromyzon marinus, Callorhinchus milii, Xenopus laevis, Mus musculus, and Anolis carolinensis). Homologs were selected according to the standard parameters used in PhylomeDB and MetaPhORs (E-value threshold < 10^–5 and a continuous overlap of 50% over the query sequence), as stated in Altenhoff et al. (2016). The number of hits included in a tree was limited to the closest 150 homologs per gene against any of the analyzed proteomes. Subsequently, alignments of homologous sequences were built in forward and reverse sense with three sequence alignment programs: MUSCLE (Edgar, 2004), MAFFT (Katoh et al., 2005), and KALIGN (Lassmann and Sonnhammer, 2005). The six resulting alignments were then combined in a consistency framework as implemented in M-COFFEE (Wallace et al., 2006), and the resulting alignment was trimmed with trimAl (consistency cut-off of 0.16667 and -gt > 0.1) (Capella-Gutiérrez et al., 2009). Multiple trees were then built, and the programing toolkit ETE (Huerta-Cepas et al., 2010) was used for each tree to understand duplication and speciation relationships by means of a 0-score species overlap phylogenetic approach. When more than one homologous was present in the same species (i.e., paralogous genes), a gene duplication was inferred and dated at the most recent common ancestor of these paralogous genes through the different lineages (Huerta-Cepas and Gabaldón, 2011). All information about orthology and paralogy relationships is available in phylomeDB (Huerta-Cepas et al., 2014). PhyML v3 (Guindon et al., 2010) was used to create a maximum likelihood tree with one-to-one orthologous in each of the selected species. Branch support was analyzed using a parametric approximate likelihood ratio test (aLRT) based on a chi-square distribution with three rates categories in all the cases. A super-tree from all single gene trees in the gilthead sea bream phylome was also reconstructed using a gene tree parsimony strategy as implemented in duptree (Wehe et al., 2008).

A functional analysis of gene ontology (GO) terms and metabolic pathways was performed over the protein coding genes (PCGs) model. Cellular Component, Molecular Function, and Biological Process GO terms were obtained from this functional analysis and a threshold of 50 counts was used to achieve the most representative GO terms for each category. Fisher test-based functional enrichment of biological process-associated GO terms was computed by analyzing the fraction of the model corresponding to chimeric/composite genes. Enrichment analysis derived from phylogenomics was also performed using FatiGO (Al-Shahrour et al., 2007) by comparing ontology annotations of the proteins involved in duplication against all the others encoded in the genome.

RNA-seq sequenced reads were processed to generate a gene expression Atlas across tissues. Briefly, reads were independently mapped against the reference transcriptome created from the set of ab initio predictions using Bowtie2. As a highly conservative procedure, only predictions with > 50% homology overlapping and ≥ 5 counts were accepted and included as reliable features. Corset v1.07 (Davidson and Oshlack, 2014) was used to quantify genes in each sample separately. Expression values were calculated in reads per kilobase per million mapped reads (RPKM) (Mortazavi et al., 2008).

To retrieve and annotate duplication events, we considered both the species-specific set of homologous genes from the phylogenomics analysis as well as ab initio predictions supported by RNA-seq transcripts. The identity of paralogs was supported by phylogenomic evidence and a same molecular description based on sequence similarity. To consider a set of paralogs as tissue-specific, all the copies must show detectable expression in more than one tissue. On the other hand, a given set of paralogs was considered as tissue-exclusive when showing expression in only one of the analyzed tissues. Student’s t-test and a one-way ANOVA (P < 0.05) were used to detect the differential expression between specialized gilthead sea bream paralogs of skeletal muscle, liver, gills and spleen. Correction by false discovery rate (FDR) (α = 0.05) was applied for all the paralog sets. This statistical analysis was not applied to intestine samples because the expression analysis was conducted with pooled instead of individual samples.

The existence of Atlas of expression in humans and other higher vertebrates^{5 ,6} was exploited to retrieve and compare the enrichment of tissue-exclusive paralogs. Accordingly, tissue-exclusive genes with non-redundant descriptions (initially assessed by RNA-seq) were categorized as follows: (1) enriched genes in the same tissue in other animal models, (2) enriched genes in the same tissue and in other tissues present in the analysis, (3) genes expressed in almost all the analyzed tissues, and (4) unclassified genes.

Duplicated genes from the analyzed tissues, covering a wide range of expression level among copies, were chosen for real-time qPCR validation: cav3, myod1, and myod2 (skeletal muscle); slc6a19 and aoc1 (intestine); upp2 and prom1 (liver); lmo1 and yjefn3 (gills); gp2 and hbb2 (spleen). GenBank accession numbers of the aforesaid duplicated transcripts are MN131091–MN131112. To complete the range of expression, cdh15 (skeletal muscle), cldn15 (intestine), clec10a (liver), sox3 (gills), and lgals1 (spleen) were included in the qPCR. The validation was performed on the same RNA individual samples used for RNA-seq. Primer design (Supplementary Table S1), reverse transcription, qPCR optimization and reactions were performed as previously detailed (Benedito-Palos et al., 2016). Specificity of reactions was verified by melting curves analyses and expression data were normalized to β-actin using the delta delta Ct method (Livak and Schmittgen, 2001). Pearson correlation coefficients were calculated in order to compare gene expression values for RNA-seq samples and qPCR expression data.

Gilthead sea bream genome was assembled using a hybrid strategy involving Illumina NextSeq500 and PacBio RS II as sequencing platforms. An overview of the main stages and achievements of the project is shown in Figure 1. Data obtained from the two PE and two MP Illumina libraries reached ∼94.8 Gb and ∼11.7 Gb, respectively (see Supplementary Table S2). PE read assembly yielded 51,918 contigs with an N50 of 50.2 kb and an L50 of 6,823 contigs. The initial assembly was further improved by means of scaffolding with MP and SMRT reads followed by gap filling. This procedure resulted in 5,039 scaffolds (>750 bp length) with an N50 scaffold length of 1.07 Mb and an L50 scaffold count of 227. At this end, the percentage of assembly in scaffolded contigs was 99.2% with a mean scaffold size of 247.38 kb and an average GC content of 39.82%. For more details in assembly metrics see Supplementary Table S3.

K-mer analysis using PE reads (Supplementary Figure S1A) showed 63-mer read length frequency with an estimated genome size of ∼1.59 Gb (main peak), including 543 Mb of repeated k-mers (repeat peak). The total scaffold length was ∼1.24 Gb, which represents 78% of the estimated total genome size. According to this, the average assembly coverage was 67.8x, and 90% of the total assembled genome was included in the largest 1,613 scaffolds (Supplementary Figure S1B).

Super-scaffolding assembly was performed using 7,700 CNEs derived from the genetic linkage map of the first gilthead sea bream genome release (Pauletto et al., 2018). These CNEs, associated to unique positions within 932 scaffolds, served for ordering and orienting 57.8% of the scaffold assembly length (∼732 Mb) in 24 super-scaffolds (Supplementary Figure S2). The resulting virtual gilthead sea bream karyotype can be viewed at www.nutrigroup-iats.org/seabreamdb.

A first ab initio prediction of PCG was carried out using AUGUSTUS v3.3 (Stanke et al., 2008). To support the establishment of the PCG model, eight RNA-seq libraries from this study (6 skeletal muscle, 2 intestine) in combination with additional libraries from liver (4), spleen (3) and gills (3) (retrieved from SRA archive) were processed to generate an Atlas of gene expression across tissues (see section “Data Availability Statement”). The sequenced reads were mapped against ab initio predictions, and 55,423 PCG were inferred based on RNA-seq transcriptome analysis and homology against SWISSPROT, NR, or the IATS-CSIC gilthead sea bream transcriptome database (Calduch-Giner et al., 2013). This procedure generated a total of 21,275 unique gene descriptions with 9,250 single-copy genes. Up to 90% of unique gene descriptions are comprised in the 1,613 largest scaffolds (Figure 2A). The average gene length is 10,134 bp with exon and intron mean sizes of 184 and 1,751 bp, respectively. This yields an average protein length of 375 amino acids. For super-scaffolded genes (i.e., genes present within the 24 super scaffolds), the number of non-redundant protein descriptions decreases to 16,046 with an average gene size of 11,756 bp (Figure 2B). Dedupe redundancy analysis performed over the transcript set retrieved from RNA-seq revealed a total of 559 genes, which represents a small fraction (1.01%) of segmental gene duplications (Supplementary Table S4). The number of containments (i.e., shorter overlapping contained sequences that could potentially result in a functional truncated protein or a false positive) at 98, 95, and 90% of identity threshold was also very low (3.31, 5.05, and 6.83%, respectively).

At the scaffold level, the gilthead sea bream mobilome accounts for the 75% of the full genome size (944 Mb). More than 60% of this mobilome (599 Mb) is constituted by introns, whereas the rest of MGEs are widely spanned throughout the assembly (Supplementary Table S5). The predicted low complexity repeats (16.91%) spanned 160.5 Mb with approximately 160 Mb corresponding to 2,500 repeat families classified as de novo specific of gilthead sea bream. The remaining 0.5 Mb corresponded to known repeats (inverted and/or tandem repeats as well as satellites and microsatellites) also present in other fish genomes. Class I MGE (5.84%) comprised 27.2 Mb of LTRs retroelements (Ty3/Gypsy, BEL/Pao, Ty1/Copia, and Retroviridae-like), 27.8 Mb of non-LTR retroelements (distributed in 14 families, mainly LINEs and SINEs), and 0.2 Mb of YR-like DIRS retrotransposons. Class II MGE (10.55%) included 99.6 Mb split in 27 groups of DNA transposons (mainly hAT, Tc1/mariner, PIF/Harbinger, and PiggyBac elements). The last fraction of the mobilome corresponded to non-coding RNA (1.25%) and chimeric/composite genes (1.95%). A complete list of non-coding RNA (ncRNA) genes is shown in Supplementary Table S6, including both long (11 Mb constituted by 10 groups; mainly lincRNA, pseudogenes, and processed transcripts) and small (1 Mb split in 11 groups mainly microRNA, tRNA, and snoRNA) ncRNA. Chimeric/composite genes (i.e., those carrying exon traits constituted by MGEs) were split in 10 groups of loci: non-LTR retroelement traits (7 Mb), LTR retroelement traits (0.7 Mb), DNA transposon traits (5.8 Mb), ncRNA gene traits (0.053 Mb), repeats (0.001 Mb), viral-related traits (0.2 Mb) and YR retroelement traits (0.02 Mb), as well as clan AA peptidases (0.047 Mb), Scan/Krab genes (0.008 Mb), and unknown genes (4 Mb). For more specific details about chimeric/composite gene annotation see Supplementary Table S7. Krona representation of split sublevels of mobilome can be seen in Supplementary Figure S3.

Functional annotation of gilthead sea bream genes using GO resulted in a diverse set of functional categories allocated to 43,221 genes (Cellular Component, 41,423; Molecular Function, 38,505; Biological Process, 38,588). The top 12 categories of each ontology for non-redundant protein descriptions are shown in Figure 3A. Cellular component GO terms had the higher gene count with cytoplasm (GO:0005737; 20,689), plasma membrane (GO:0005886; 16,138) and integral to membrane (GO:0016021; 12,436) GO terms. The most abundant Molecular Function GO terms comprised metal ion binding (GO:0043167; 9,210), DNA binding (GO:0003677; 7,041) and ATP binding (GO:0005524; 6,518). The most represented biological process GO terms were transcription DNA-dependent (GO:0006351; 6,222), signal transduction (GO:0007165; 3,851) and multicellular organismal development (GO: 0007275; 2,908).

When tested for enrichment of GO terms among chimeric/composite genes, the 3,648 duplicated genes with 108 non-redundant protein annotations (Supplementary Table S8) rendered 184 enriched biological processes (corrected P-value < 0.05). These genes cover different GO terms related to immune system (26%), cell cycle (16%), translational initiation (11%), response to activity (11%), signal transduction (6%), developmental process (5%), and growth (2%) among others (Figure 3B). The relationship among functional categories is illustrated by a Venn diagram, showing 87 non-redundant gene descriptions of the main five functional categories (Figure 3C). This procedure highlighted that the high representation of immune system in chimeric/composite genes was mostly due to a wide overlapping of immune GO terms with the other enriched functional categories. Intriguingly, main intersections were found among immune system process, cell cycle and signal transduction, comprising 15 enriched GO terms and 15 unique gene descriptions, corresponding to different isoforms of protein NLRC3 and NACTH, LRR, and PYD domains-containing protein 12.

Homology relationships between genes contained in the assembled gilthead sea bream super-scaffolds and genes sequenced in other species, as well as their syntenic relationships were studied. From the 30,455 gilthead sea bream genes included in super-scaffolds, 25,806 (84.73%) had orthologs in at least one of the analyzed species, being Nile tilapia (O. niloticus, 20,561), zebra mbuna (M. zebra, 19,717), platyfish (X. maculatus, 15,093), and stickleback (G. aculeatus, 14,612) the species sharing more orthologous genes with gilthead sea bream, whereas the lowest numbers of orthologous were obtained in rainbow trout (8,866) and zebrafish (D. rerio, 4,288) (Figure 4A). Likewise, the number of syntenic blocks ranged between 483 in O. niloticus to 32 in D. rerio (Supplementary Table S9). Thus, the levels of both orthology and synteny conservation reflects phylogenetic proximity among the compared species. Also, the number of orthologous genes in syntenic blocks were maximal in O. niloticus (9,914; 30.02%), M. zebra (9,499; 34.48%) and G. aculeatus (6,866; 46.85%), whereas salmonids and cyprinids showed the lowest levels of synteny with 1,284 (O. mykiss), 1,482 (Atlantic salmon, S. salar) and 44 (D. rerio) orthologous in syntenic blocks. The intra-species synteny analysis (i.e., search of syntenic blocks in the own gilthead sea bream genome) rendered a total of 268 syntenic blocks that comprised 1,131 paralogs. This large number of syntenic genes as well as the high number of connections (homology between gilthead sea bream genes) in the Circos plot of Figure 4A are indicative of a highly duplicated genome.

To gain insights in the evolution of gilthead sea bream genome and study in more detail the origin of these high levels of genomic duplication, we inferred its phylome -i.e., the complete collection of gene evolutionary histories- across 19 fully-sequenced vertebrate species. To provide a phylogenetic context to our comparisons, we reconstructed a species tree. This was made using two complementary approaches: (1) species tree concatenation of a total of 148 genes with one-to-one orthologous in each of the included species and (2) super-tree reconstruction using 58,484 gene trees from the phylome. Both approaches resulted in the same highly supported topology (Figure 4B), which was fully consistent with the known relationships of the considered species. All trees and alignments are available to browse or download through PhylomeDB⁷ (Huerta-Cepas et al., 2014) under the phylomeDB ID 714.

From the reconstructed gilthead sea bream phylome, we inferred that 45,162 genes had duplications. The fraction of duplicated genes remained high (17,596) after the removal of gene family expansions (i.e., those resulting in five or more in-paralogs). When duplication frequencies per branch in all lineages leading to the gilthead sea bream were computed, two peaks of high duplication ratios (average duplications per gene) were inferred at earliest splits of vertebrates and at the base of teleost fish (teleost-specific genome duplication), which correspond to the known WGDs (Figure 4B; clades 8, 12). Additionally, the gilthead sea bream genome also showed a high rate of species-specific duplications (2.024 duplications per gene; 0.385 duplications per gene after removing expansions). Functional GO enrichment of these duplicated genes highlighted different biological processes, mostly related to genome transposition, immune response and response to stimulus. This referred to the following GO terms: DNA integration (GO:0015074); transposition, DNA-mediated (GO:0006313); RNA-dependent DNA biosynthetic process (GO:0006278); developmental process (GO:0032502); transposition, RNA-mediated (GO:0032197); DNA recombination (GO:0006310); immunoglobulin production (GO:0002377); detection of chemical stimulus involved in sensory perception (GO:0050907); regulation of T cell apoptotic process (GO:0070232); telomere maintenance (GO:0000723). In the case of immunoglobulin production, this stated to 24 unique gene descriptions including among others Ig heavy chain Mem5-like isoform X1, Ig heavy chain Mem5-like isoform X2, Ig kappa chain V region 3547, Ig kappa chain V region Mem5, Ig kappa chain V-II region 2S1.3, Ig kappa chain V-IV region Len, Ig lambda chain V-I region BL2, Ig lambda chain V-I region NIG-64, Ig lambda-3 chain C regions, Ig lambda-6 chain C region, Ig lambda-6 chain C region, Ig lambda-like polypeptide 1 isoforms X1, X3 and X4, Ig lambda-like polypeptide 5, pre-B lymphocyte protein 3, integral membrane protein 2A, laminin subunit alpha-2, or Ig kappa chin V19-17. Likewise, the regulation of T cell apoptotic process refers to microfibrillar-associated protein 1, tyrosine-protein kinase JAK2 and JAK3 in addition to different GTPases of IMAP family members (2, 4, 4-like, 8, 8-like). Lastly, the category detection of chemical stimulus involved a wide representation of olfactory receptors, including among others olfactory receptor 10J4-like, 11A11-like, 13C8-like, 146-like, 1M1-like, 2K2-like, 2S2-like, 4C15-like, 4K3-like, 4N5-like, 51G1-like, 5A5-like, 52D1-like, 52K1-like, 5B17-like, and 6N1-like.

Up to 70% of the pre-processed reads of the RNA-seq tissue samples were mapped in the assembled genome, yielding 55,423 genes that are reduced to 16,992 after the removal of low expressed genes, low alignments high scoring pairs (HSPs) and phylome-based paralogs. From these filtered sequences, up to 5,322 genes were recognized as ubiquitously expressed sequences in the analyzed tissues (Figure 5A). Intestine as a whole (anterior and posterior intestine segments) had the highest number of tissue-exclusive annotated genes (1,198), followed by gills (667), liver (256) and spleen (248) and skeletal white muscle (203). When unique gene descriptions were considered, the order of tissues with a tissue-exclusive number of non-redundant molecular signatures was maintained: intestine (512) > gills (379) > liver (139) > spleen (131) > skeletal muscle (123) (Figure 5B). This yielded a variable percentage of duplicated genes from 28% in the consensus gene list (1.295 out of 4.625) for all the analyzed tissues to 20–17% in muscle and intestine, 12–10% in liver and gills and 6% in spleen. Likewise, the duplication rate ranged between 1.62 from the consensus list to 1.26–1.24 in muscle and intestine, 1.16 in liver, 1.13 in gills, and 1.08 in spleen (Figure 5C). The final list of 1,284 tissue-exclusive genes (present in only one tissue) with their number of copies is shown in Supplementary Table S10.

Tissue-exclusive non-redundant paralogs of intestine, skeletal muscle, liver, spleen, and gills are listed in Supplementary Table S11. According to the gene expression pattern in humans and other higher vertebrates^5,6, most of them (65–75%) were classified as tissue- or group-enriched genes (gills paralogs are not included in the analysis due to the lack of a reference expression Atlas for fish species) (Figure 6A). This procedure yielded up to 65 tissue-exclusive paralogs (intestine, 30; skeletal muscle, 17; liver, 13; spleen, 5), showing expression changes between duplicated copies with a similar range of variation when the outliers from intestine (1) and gills (1) were not included in the analysis (Figure 6B). For some of them, including cav3, myod1 and myod2 (skeletal muscle); slc6a19 and aoc1 (intestine); upp2 and prom1 (liver); lmo1 and yjefn3 (gills); gp2 and hbb2 (spleen) the differential gene expression pattern for duplicated genes was validated by qPCR, and overall a high correlation was found for representative genes of all analyzed tissues (Supplementary Table S12).