Genome features

We applied a bacterial artificial chromosome (BAC) and whole-genome shotgun (WGS) hierarchical assembly strategy for the L. crocea genome to overcome the high levels of genome heterozygosity (Table 1 and S1–S2 Fig). The 42,528 BACs were sequenced by the HiSeq 2000 platform and each BAC was assembled by SOAPdenovo [18] (S1 Table). The total length of all combined BACs was 3,006 megabases (Mb), which corresponded to approximately 4.3-fold genome coverage (S2–S3 Tables). All BAC assemblies were then merged into super-contigs and oriented to super-scaffolds with large mate-paired libraries (2–40 kb). Gap filling was made with reads from short insert-sized libraries (170–500 bp) (S3–S4 Tables). In total, we sequenced 563-fold coverage bases of the estimated 691 Mb genome size. The final assembly was 679 Mb, with a contig N50 of 63.11 kb and a scaffold N50 of 1.03 Mb (Table 1). The 672 longest scaffolds (11.2% of all scaffolds) covered more than 90% of the assembly (S5 Table). To assess the completeness of the L. crocea assembly, 52-fold coverage paired-end high-quality reads were aligned against the assembly (S3 Fig). More than 95.63% of the generated reads could be mapped to the assembly by Burrows-Wheeler Aligner and 98.99% of assembled sequences could be covered by at least 5 reads. Furthermore, the integrity of the assembly was validated by the successful mapping of 95.80% of 18,184 transcripts (greater than 1000 bp) with an identity cutoff of 80% (S6 Table). These results indicate that the genome assembly of L. crocea has high coverage and is of high quality (S7 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of the Larimichthys crocea genome

https://doi.org/10.1371/journal.pgen.1005118.t001

The repetitive elements comprise 18.1% of the L.crocea genome (S8 Table), which is a relatively low percentage when compared with other fish species, such as Danio rerio (52.2%), Gadus morhua (25.4%), and Gasterosteus aculeatus (25.2%). This suggests that L. crocea may have a more compact genome (S9–S10 Tables). We identified 25,401 protein-coding genes based on ab initio gene prediction and evidence-based searches from the reference proteomes of six other teleost fish and humans (S4 Fig and S11 Table), in which 24,941 genes (98.20% of the whole gene set) were supported by homology or RNAseq evidence (S5 Fig). Over 97.35% of the inferred proteins matched entries in the InterPro, SWISS-PROT, KEGG or TrEMBL database (S12 Table).

Phylogenetic relationships and genomic comparison

The phylogenetic relationships of L. crocea to seven other sequenced teleost species were estimated based on 2,257 one-to-one high-quality orthologues, using the maximum likelihood method. According to the phylogeny and the fossil record of teleosts, we dated the divergence of L. crocea from the other teleost species to approximately 64.7 million years ago (Fig. 1A). We also detected 19,283 orthologous gene families (S3 Table), of which 14,698 families were found in L. crocea. The gene components of L. crocea were similar to those of D. rerio (Fig. 1B). The gene contents in four representative teleost species and L. crocea genomes were also analysed, and 11,205 (76.23%) gene families were found to be shared by five teleosts (Fig. 1C). We confirmed that the one-to-one orthologous genes of G. aculeatus and L. crocea have higher sequence identities from the distribution of the percent identity of proteins (Fig. 1D), which indicates that Sciaenidae has a closer affinity to Gasterosteiformes and coincides with our genome-level phylogeny position.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Phylogenetic tree of and orthologous genes in L. crocea and other vertebrates.

(A) The phylogenetic tree was constructed from 2,257 single-copy genes with 3.18 M reliable sites by maximum likelihood methods. The red points on six of the internal nodes indicate fossil calibration times in the analysis. Blue numbers indicate the divergence time (Myr, million years ago), and the green and red numbers represent the expanded and extracted gene families, respectively, in L. crocea. (B) The different types of orthologous relationships are shown. “1:1:1” = universal single-copy genes; “N:N:N” = orthologues exist in all genomes; “Fish” = fish-specific genes; “SD” = genes that have undergone species-specific duplication; “Homology” = genes with an e-value less than 1e-5 by BLAST but do not cluster to a gene family; “ND” = species-specific genes; and “Others” = orthologues that do not fit into the other categories. (C) The shared and unique gene families in five teleost fish are shown in the Venn diagram. (D) Distribution of the identity values of orthologous genes is compared among L. crocea and other teleosts.

https://doi.org/10.1371/journal.pgen.1005118.g001

Furthermore, 121 significantly expanded and 27 contracted gene families (P < 0.01) were identified by comparing the family size of L. crocea with that of the other vertebrates used in the phylogenetic analysis (S14–S15 Tables). Based on the ratios of the number of nonsynonymous substitutions per nonsynonymous site (K_a) to the number of synonymous substitutions per synonymous site (K_s; K_a/K_s ratios) in a branch-site model of PAML [19], 92 genes were found to be positively selected in L. crocea compared with their orthologues in the other six teleost species (P < 0.001, S16 Table).

Genetic features of the L. crocea

L. crocea is a migratory fish with good photosensitivity, olfactory detection, and sound perception, and it contains high levels of selenium [8]. Our genomic analyses provide genetic basis for these physiological characteristics. Several crystallin genes (crygm2b, cryba1, and crybb3), which encode proteins that maintain the transparency and refractive index of the lens [4], were significantly expanded in the genome of L. crocea relative to those of other sequenced teleosts (S17 Table). Phylogenetic analysis showed that the crystallin genes from L. crocea cluster together, indicating that these genes were specifically duplicated in L. crocea lineage (S6 Fig). The specific expansion of these crystallin genes may be helpful for improving photosensitivity by increasing lens transparency, thereby enabling the fish to easily find food and avoid predation underwater.

We also identified 112 olfactory receptor (OR)-like genes from the L. crocea genome (S18 Table and S7 Fig.), and almost all of them (111) have been reported to be expressed in the olfactory epithelial tissues of L. crocea [14]. The majority of these genes (66) were classified into the “delta” group, which is important for the perception of water-borne odorants [20]. L. crocea also possessed the highest number of genes that were classified into the “eta” group (30, P < 0.001), and these genes may contribute to the olfactory detection abilities, which could be useful for feeding and migration [21].

L. crocea is named for its ability to generate strong repetitive drumming sounds, especially during reproduction [8]. For good communication, fish have developed high sensitivities to environmental sound. Three important auditory genes, otoferlin (OTOF), claudinj, and otolin 1 (OTOL1), were significantly expanded in the L. crocea genome (P < 0.01, S19 Table). These expansions may contribute to the detection of sound signaling during communication, and thus to reproduction and survival [22].

Selenium is highly enriched in L. crocea [8], and it is mainly present as selenoproteins. We used the SelGenAmic-based selenoprotein prediction method [23] to analyse the L. crocea genome and identified 40 selenoprotein genes, which is the highest number among all examined vertebrates (S20 Table). Interestingly, five copies of MsrB1, which encodes methionine sulfoxide reductase, were found in L. crocea (MsrB1a, MsrB1b, MsrB1c, MsrB1d, and MsrB1e), whereas only two copies (MsrB1a and MsrB1b) were found in other fish, thus suggesting its broader specificity to reduce all possible substrates [24].

Characterization of the L. crocea immune system

Approximately 2,528 immune-relevant genes were annotated in the L. crocea genome, including 819 innate immune-relevant genes and 1,709 adaptive immune-relevant genes (S21 Table). Strikingly, L. crocea was found to have not only a relatively complete innate immune system, but also a well-established adaptive immune system, because the majority of the CD4⁺ T-helper type 1 (Th1), CD4⁺ T-helper type 2 (Th2) and CD8⁺ T cell-related genes were found (Fig. 2A). This suggests that the CD4⁺ Th1, CD4⁺ Th2 and CD8⁺ T cell-mediated adaptive immunity is well conserved in L. crocea. Moreover, the genes related to Th17 cell- and γδ-T cell-mediated mucosal immune responses were also conserved in L. crocea, suggesting that a well-developed mucosal immunity exists in this species as well (Fig. 2A). We detected gene expansions in several of these immune-relevant genes, including those encoding lectin receptors (CLEC17A), a classical complement component (C1q), apoptosis regulator (BAX), and immunoglobulins (IgHV) (P < 0.01, S22 Table). Expansions were also observed in the genes encoding four key proteins for mammalian antiviral immunity: tripartite motif containing 25 (TRIM25), cyclic GMP-AMP synthase (cGAS), DDX41, and NOD-like receptor family CARD domain containing 3 (NLRC3) (Fig. 2B). However, retinoic acid-inducible gene-1 (RIG-I), which initiates antiviral signaling pathway by recognizing cytosolic virus-derived RNA in mammals, was not found in the L. crocea genome and transcriptome [13,16]. The teleost RIG-I has been identified only in limited fish species, such as cyprinids and salmonids, and its absence suggests that it may have been lost from particular fish genomes [25]. In mammals, laboratory of genetics and physiology 2 (LGP2) can serve as a suppressor to block RIG-I- and melanoma differentiation-associated protein 5 (MDA5)-elicited signaling [26]. However, LGP2 in fish can bind poly (I:C) to trigger interferon (IFN) production [27], whose functional performance was similar to that done by RIG-I, thereby, LGP2 in fish may act as a substitute for RIG-I (Fig. 2B). The expanded TRIM25 (54 copies, S8 Fig) may trigger the ubiquitination of IFN-β promoter stimulator-1 (IPS-1), thus allowing for IFN regulatory factor 3 (IRF3) phosphorylation and antiviral signaling initiation [28]. DDX41 and cGAS encode intracellular DNA sensors that activate stimulator of interferon genes (STING) and TANK-binding kinase 1 (TBK1) to induce production of type I IFNs [29,30]. Fish type I IFNs can be divided into two major groups, namely group I and II type I IFNs, with two or four cysteines forming one or two pairs of intramolecular disulphide bonds, respectively [31]. In the L. crocea genome, two type-I IFN genes, tentatively named IFN-d and IFN-h, were detected, belonging to group I type- I IFNs (Fig. 2B). L. crocea contained 43 copies of NLRC3 (S9 Fig), which encodes regulators that prevent type I IFN overproduction [32]. The expansions of these virus-response genes suggest their enhanced roles in innate antiviral immunity, which may explain why L. crocea is less susceptible to viral infection.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Characterisation of the T-cell lineages in L. crocea adaptive immunity and the expanded genes in antiviral immunity.

(A) A schematic diagram summarising genes related to different T-cell lineages in L. crocea is shown. The inducible factors, the main regulatory transcriptional factors, and the immune effectors of T cells are present in green, blue, and orange backgrounds, respectively. The genes that have been annotated by genome survey are shown in black, and the unannotated genes are shown in red. IL-2R represents all three subunit genes of L. crocea IL-2 receptor, IL-2RB, IL-2RG and IL-2RA/ IL-15RA genes. The majority of the CD4⁺ T-helper type 1 (Th1), CD4⁺ T-helper type 2 (Th2) and CD8⁺ T cell-related genes have been found in the L. crocea genome. The genes related to Th17 cell- and γδ-T cell-mediated mucosal immunity are also conserved in L. crocea. (B) Several key genes are expanded in the antiviral immunity pathways in L. crocea. The genes that have been identified in the L. crocea genome are shown in orange boxes, and the lost gene (RIG-I) is shown in the grey box. LGP2 is able to bind to double-stranded RNA (dsRNA) to trigger interferon production, but the adaptor molecule of LGP2 is still unknown in fish. The red boxes indicate gene families (TRIM25, cGAS, DDX41, and NLRC3) that are expanded in L. crocea. The arrow represents induction, and the interrupted line represents inhibition.

https://doi.org/10.1371/journal.pgen.1005118.g002

Stress response under hypoxia

The brain allows rapid and coordinated responses to the environmental stress by driving the secretion of hormones. Therefore, we studied the response of the L. crocea brain to hypoxia. We sequenced seven transcriptomes of the brains at different times of hypoxia exposure and found that 8,402 genes were differentially expressed at one or more time points (fold change ≥ 2, false discovery rate [FDR] ≤ 0.001; S10 Fig). Hypoxia stress induced a response with the largest number of genes (4,535 genes) at 6 h (S11 Fig), indicating that genes with regulated expression at 6 h may be critical for the response.

Among the numerous incidents happened during the early hypoxia stress, activation of the central neuroimmune system has been shown the most critical event, in which brain neuropeptides, endocrine hormones, and inflammatory cytokines closely participate to generate protective effects [33–35]. However, the precise regulatory networks among these factors have not yet been fully delineated. Our transcriptome analyses show that the key hypothalamic-pituitary-adrenal (HPA) axis-relevant genes (corticotropin-releasing factor [CRF], CRF receptor 1 [CRFR1], pro-opiomelanocortin [POMC], and CRF-binding protein [CRFBP]) in the L. crocea brain displayed a down-up-down-up (W-type) dynamic expression pattern under hypoxia stress (S12 Fig). In contrast, the endothelin-1 (ET-1) and adrenomedullin (ADM) genes showed an up-down-up-down (M-type) dynamic expression pattern, and the time of inflexion point corresponded with that of CRF, CRFR1, POMC, and CRFBP. The HPA axis can strictly control the production of glucocorticoids [36,37], and glucocorticoids are suppressors of ET-1 and ADM, both of which are involved in cerebral inflammation in mammals [38,39]. We therefore suggest that a feedback regulatory pathway may exist between the HPA axis and ET-1/ADM under hypoxia stimulation. As a support, the expression of the inflammatory cytokine genes (IL-6/TNF-α) also showed the M-type pattern and was consistent with that of ET-1/ADM (S12 Fig). Expression changes of selected genes, including CRF, CRFR1, POMC, ET-1, ADM and IL-6, during hypoxia stress, were further confirmed by real-time PCR and found to be well corresponding to their expression changes in the transcriptome analyses (S13 Fig). These coordinated and fluctuating expression patterns indicate that hypoxia may inhibit the HPA axis and induce the expression of ET-1/ADM, and the latter then promote the expression of IL-6/TNF-α and form a positive feedback loop with them (Fig. 3). On the other hand, the ET-1/ADM and IL-6/TNF-α may in turn activate the HPA axis, and the latter subsequently induces glucocorticoids and generates a negative feedback to inhibit ET-1/ADM and IL-6/TNF-α expression to reduce the over-inflammatory responses in brain, the latter of which has been reported in mammals [39–42]. Therefore, our results may outline a novel HPA axis-ET-1/ADM-IL-6/TNF-α feedback regulatory loop in the neuro-endocrine-immune network during hypoxia responses (Fig. 3 and S12–S13 Fig). Besides, we also found that both SOCS-1 and SOCS-3 in the L. crocea brain display opposite expression patterns against IL-6 and TNF-α (S12–S13 Fig). Thus, SOCS-1 and SOCS-3 may have complementary roles in down-regulating IL-6 and TNF-α, and both IL-6 and TNF-α have reciprocal functions to induce SOCS-1 and SOCS-3 expression (Fig. 3). These results suggest that a SOCS-1/3-dependent feedback regulation may exist in the process against hypoxia-induced cerebral inflammation in L. crocea.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Hypoxia stress exerts responses involving the HPA and HPT axes.

Under hypoxia, the potential neuro-endocrine-immune/metabolism networks contribute to the regulation of moderate inflammation and the maintenance of energy balance. Hypoxia may inhibit the hypothalamic-pituitary-adrenal (HPA) axis and induce the expression of ET-1/ADM, and the latter then promote the expression of IL-6/TNF-α and form a positive feedback loop with them. On the other hand, the ET-1/ADM and IL-6/TNF-α may in turn activate the HPA axis, and the latter subsequently induces glucocorticoids and generates a negative feedback to inhibit ET-1/ADM and IL-6/TNF-α expression to reduce the inflammatory responses in brain. Therefore a HPA axis-ET-1/ADM-IL-6/TNF-α feedback regulatory loop is outlined here. Besides, both SOCS-1 and SOCS-3 may have complementary roles in down-regulating IL-6 and TNF-α, and both IL-6 and TNF-α have reciprocal functions to induce SOCS-1 and SOCS-3 expression, thus constituting a SOCS-1/3-dependent feedback regulation to modify cerebral inflammation. The hypothalamic-pituitary-thyroid (HPT) axis was inhibited in L. crocea brains during the early period of hypoxia, thus leading to a decrease in thyroid hormone (T3 and T4) production. Decrease of HPT axis-thyroid hormones subsequently may inhibit protein synthesis by the PI3K-Akt-mTOR-S6K signaling pathway, which contributes to the reduction in energy consumption during hypoxia stress. Meanwhile, down-regulation of HPT axis-thyroid hormones also represses the tricarboxylic acid (TCA) cycle and accelerates the anaerobic glycolytic pathway, which facilitates O₂-independent ATP production under hypoxia. Therefore the HPT axis-mediated effects may play roles in response to hypoxia by reorganizing energy consumption and generation. Genes related to the neuro-endocrine system (orange), immunity (red), and metabolic system and protein synthesis (blue) are indicated. The outer border indicates the brain of L. crocea. The arrow represents promotion, and the interrupted line represents inhibition. Solid lines indicate direct relationships between genes. Dashed lines indicate that more than one step is involved in the process.

https://doi.org/10.1371/journal.pgen.1005118.g003

Energy maintenance and homeostasis under hypoxia is another critical event in brain. However, the mechanisms underlying this issue remain limited. Here, we found that the major hypothalamic-pituitary-thyroid (HPT) axis-related genes (thyrotropin-releasing hormone [TRH], TRH receptor [TRHR], thyroid-stimulating hormone β [TSHβ], and thyroid hormone receptor α [TRα]) were significantly down-regulated in the L. crocea brain at 1 h to 6 h under hypoxia (S14 Fig). Down-regulation of these gene expressions during hypoxia was also confirmed by real-time PCR (S15 Fig). These findings therefore indicate that the HPT axis may be inhibited during the early period of hypoxia, which agreed with the previous observations that hypoxia could influence the HPT axis [43].

HPT axis can regulate protein synthesis and glucose metabolism by production of thyroid hormones [44,45]. Inhibition of the HPT axis leads to a decrease in the production of thyroid hormones. Thyroid hormones can regulate ribosomal biogenesis and protein translation by the PI3K-Akt-mTOR-S6K signaling pathway [44]. In this study, the mRNA levels of PI3K, S6K, and most of the components of the protein translation machinery, including the ribosomal proteins and eukaryotic translation initiation factors (eIF-1, -2, -3, -5 and -6), were all down-regulated under hypoxia (S16 Fig). This suggests that the HPT axis may inhibit protein synthesis under hypoxia by decreasing the production of thyroid hormones (Fig. 3), which is beneficial for saving energy during hypoxia stress. Thyroid hormones can also accelerate the oxidative metabolism of glucose and inhibit the glycolytic anaerobic pathway [46]. Interestingly, genes involved in the tricarboxylic acid (TCA) cycle (pyruvate dehydrogenase complex [PDC-E1], succinyl-CoA synthetase [SCS], and fumarate hydratase [FH]) were down-regulated 12 h later under hypoxia, whereas glycolysis-related genes, such as pyruvate kinase (PKM), glyceraldehyde 3-phosphate dehydrogenase (GAPDH), GPI, and aldolase A (ALDOA), were greatly increased at 1 h respectively (S14–S15 Fig). The down-regulation of HPT axis-thyroid hormones may inhibit the TCA cycle and accelerate the anaerobic glycolytic pathway in the brain during hypoxia exposure (Fig. 3). The repression of the TCA cycle and the strong induction of the anaerobic glycolytic pathway resulted in a physiological shift from aerobic to anaerobic metabolism, where fish utilise O₂-independent mechanisms to produce adenosine triphosphate (ATP). However, the mRNA levels of hypoxia-inducible factor (HIF)-1α, which are significantly up-regulated under hypoxia in mammals [47,48], were not significantly changed in the L. crocea brain (S14 Fig). It is possible that the HIF-1α-mediated mechanism may not be essential for the hypoxia response in the L. crocea brain during the early period of hypoxia. These results suggest that the HPT axis-mediated effects may play major roles in response to hypoxia by reorganizing energy consumption and energy generation.

Mucus components and function

The skin mucus is considered as the first defensive barrier between fish and its aquatic environment, and it plays a role in a number of functions, including locomotion, antioxidant responses, respiration, disease resistance, communication, ionic and osmotic regulation [49]. However, the exact mechanisms underlying these functions remain unknown. Mucus is composed mainly of the gel-forming macromolecule mucins and inorganic salts suspended in water [50]. We identified 159 genes that are implicated in mucin biosynthesis and mucus production in the L. crocea genome (S23 Table), based on previous studies in mammals [51]. This indicates that the mucin synthetic pathway is conserved between fish and mammals. Among these gene families, GALNT, which encodes N-acetylgalactosaminyl transferases [52], was significantly expanded in L. crocea (27 copies versus 15–20 copies in other fish; S17 Fig). Syntaxin-11 was also expanded. Additionally, genes encoding syntaxin-binding protein 1 and syntaxin-binding protein 5, which are related to mucus secretion, were positively selected in the L. crocea genome (S16 Table). The expansion and positive selection of these genes may explain why the L. crocea secretes more mucus than other fish under stress.

We identified 22,054 peptides belonging to 3,209 genes in the L. crocea skin mucus proteome, and this accounted for more than 12% of the protein-coding genes in the genome (S24 Table). The complexity of the L. crocea mucus presumably relates to the multitude of its biological functions that allow the fish to survive and adapt to environmental changes. The over-represented functional categories were oxidoreductase activity (GO:0016491, P = 1.58×10^-35, 223 proteins), peroxidase activity (GO:0004601, P = 0.0075, 9 proteins), oxygen binding (GO:0019825, P = 0.0011, 8 proteins), and ion binding (GO:0043167, P = 2.21×10^-6, 347 proteins) (Fig. 4A and S18 Fig.). Two hundred and thirty-two antioxidant proteins that were related to oxidoreductase activity and peroxidase activity were highly enriched in the L. crocea mucus, and they included peroxiredoxins, glutathione peroxidase, and thioredoxin (S25 Table). These proteins intercept and degrade environmental peroxyl and hydroxyl radicals from aqueous environments [53]. Therefore, the presence of high-abundance antioxidant proteins in the skin mucus may have the potential to protect fish against air exposure-induced oxidative damage (Fig. 4B). Eight proteins related to oxygen transport, including hemoglobin subunits α1, αA, αD, β, and β1, and cytoglobin-1, were identified in the L. crocea skin mucus (S26 Table). The abundant expression of hemoglobin may contribute to the binding and holding of oxygen for respiration. Various immune molecules that provide immediate protection to fish from potential pathogens, such as lectins, lysozymes, C-reactive proteins, complement components, immunoglobulins, and chemokines, were also found in the L. crocea skin mucus (S27 Table). To date, the mechanisms of osmotic and ionic regulation of the skin mucus have not been confirmed [49]. In this study, a large number of ion-binding proteins were identified in the L. crocea mucus (S28 Table). These proteins and the layer of mucus may have a role in limiting the diffusion of ions on the surface of the fish (Fig. 4B). However, a substantial proportion of the proteins, which are highly present in the skin mucus of fish under air exposure, play an unknown role in the mucus response.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Skin mucus proteins are overexpressed in air-exposed L. crocea.

(A) The distribution of mucus proteins in the molecular function class of Gene Ontology is shown. The over-represented functional categories are indicated in the pie chart, of which oxidoreductase activity-, peroxidase activity-, oxygen binding-, and ion binding-related proteins are enriched. (B) A representation of the functional mechanisms of the mucus barrier is shown. The continuously replenished thick mucus layer can retain a large number of antioxidant, immune, oxygen-binding, and ion-binding molecules, which are involved in antioxidant functions, immune defence, oxygen transport, and osmotic and ionic regulation, respectively. Antioxidant, immune, oxygen-binding, and ion-binding molecules are indicated in red, green, blue, and gray, respectively.

https://doi.org/10.1371/journal.pgen.1005118.g004

Ethics statement

The studies were carried out in strict accordance with the Regulations of the Administration of Affairs Concerning Experimental Animals established by the Fujian Provincial Department of Science and Technology. Animal experiments were approved by the Animal Care and Use Committee of the Third Institute of Oceanography, State Oceanic Administration. All surgery was performed under Tricaine-S anesthesia, and all efforts were made to minimize suffering.

Genome assembly

The wild L. crocea individuals were collected from the Sanduao sea area in Ningde, Fujian, China. Genomic DNA was isolated from the blood of a female fish by using standard molecular biology techniques for BAC library construction and sequencing by the HiSeq 2000 Sequencing System in BGI (Beijing Genomics Institute, Shenzhen, China). Subsequently, low-quality and duplicated reads were filtered out, and sequencing errors were removed. The BACs of L. crocea were assembled by using SOAPdenovo2 [60] (http://soap.genomics.org.cn) with k-mers that ranged from 25 to 63 in size. Then, we selected the assembly with the longest scaffold N50 for gap filling. The BACs were merged together based on the overlap found by BLAT, using the custom script: Rabbit (ftp://ftp.genomics.org.cn/pub/Plutellaxylostella/Rabbit_linux-2.6.18-194.blc.tar.gz). The redundant sequences that were produced by high polymorphisms were removed by sequence depth and shared k-mer percentage. Assembly was performed by scaffolding with mate-paired libraries (2–40 kb) using SSPACE v2 [61], and gap filling was made by Gapcloser (http://sourceforge.net/projects/soapdenovo2/files/GapCloser/) with small-insert libraries (170–500 bp).

The whole-genome sequences of L. crocea have been deposited in the DNA Data Bank of Japan (DDBJ), the European Molecular Biology Laboratory (EMBL) nucleotide sequencing database and GenBank under the accession JRPU00000000. The version described in this paper is the first version JRPU01000000. All short-read data of WGS and BAC have been deposited in the Short Read Archive (SRA) under accession SRA159210 and SRA159209 respectively.

Genome annotation

For the annotation of repetitive elements, we used a combination of homology-based and ab initio predictions. RepeatMasker [62] and Protein-based RepeatMasking [62] were used to search Repbase, which contains a vast amount of known transcriptional elements at the DNA and protein levels. During the process of ab initio prediction, RepeatScout [63] was used to build the ab initio repeat library based on k-mer, using the fit-preferred alignment score on the L. crocea genome. Contamination and multi-copy genes in the library were filtered out before the RepeatScout library was used to find homologs in the genome and to categorise the found repeats by RepeatMasker [62].

Gene models were integrated based on ab initio predictions, homologue prediction, and transcription evidence.

Homology-based prediction

The protein sequences of seven species (Danio rerio, Gasterosteus aculeatus, Oreochromis niloticus, Oryzias latipes, Takifugu rubripes, Tetraodon nigroviridis, and Homo sapiens) were aligned to the L. crocea assembly using BLAST (E-value ≤ 1e-5), and the matches with length coverage >30% of the homologous proteins were considered as gene-model candidates. The corresponding homologous genome sequences were then aligned against the matching proteins by using Genewise [64] to improve gene models.

Ab initio prediction

Augustus [65], SNAP [66], and GENESCAN [67] were used for the ab initio predictions of gene structures on the repeat-masked assembly.

Transcriptome-based prediction

RNAseq reads from the transcriptomes of the mixed tissues of a female and a male (eleven tissues each) were aligned to the genome assembly by Tophat [68], which can identify splice junctions between exons. Cufflinks [69] was used to obtain transcript structures.

Homology-based, ab initio derived and transcript gene sets were integrated to form a comprehensive and non-redundant gene set. The overlap length of each gene was verified by different methods, and genes showing 50% overlap by at least one method were selected. To eliminate false positives (genes only supported by ab initio methods), novel genes with the reads per kilobase of gene model per million mapped reads (RPKM) ≤ 1 were removed. A method based on gene synteny and FGENESH program (http://www.softberry.com/berry.phtml?topic=fgenesh&group=programs&subgroup=gfind&advanced=on) was used for the prediction of interleukin genes.

Evolutionary and comparative analyses

To detect variations in the L. crocea genome, we chose nine species (Larimichthys crocea, Gasterosteus aculeatus, Takifugu rubripes, Tetraodon nigroviridis, Oryzias latipes, Gadus morhua, Danio rerio, Gallus gallus, and Homo sapiens). Proteins that were greater than 50 amino acids in size were aligned by BLAST (-p blastp-e 1e-7), and Treefam [70] was used to construct gene families for comparison.

The 2,257 single-copy genes from the gene family analysis were aligned using MUSCLE [71], and alignments were concatenated as a single data set. To reduce the error topology of phylogeny by alignment inaccuracies, we used Gblock [72] (codon model) to remove unreliably aligned sites and gaps in the alignments. The phylogenetic tree and divergence time were calculated using the PAML 3.0 package [19].

Gene family expansion and contraction analyses were performed by cafe [73]. For optical, olfactory receptor, and auditory system-related genes, we downloaded the genes from Swissprot or Genebank and predicted their candidates using BLAST and Genewise to determinate copy numbers. Pseudogenes produced by frame shift were removed. Phylogenetic analysis of the expanded gene families was based on maximum likelihood methods by PAML 3.0 [19], and the phylogenetic tree was represented by EvolView [74].

Amino acid sequences from six representative teleosts (Larimichthys crocea, Gasterosteus aculeatus, Danio rerio, Oryzias latipes, Takifugu rubripes, and Tetraodon nigroviridis) were aligned by BLAST (-p blastp-e 1e-5-m 8), and reciprocal-best-BLAST-hit methods were used to define orthologous genes in six teleost fish. Because alignment errors are an important concern in molecular data analysis, we made alignments of codon sequences, which are nucleotide sequences that code for proteins, using the PRANK aligner [75]. Positive selection was inferred, based on the branch-site K_a/K_s test by codeml in the PAML 3.0 package [19].

Transcriptome under hypoxia

L. crocea (90–100 g) individuals were purchased from a mariculture farm in Ningde, Fujian, China. The fish were maintained at 25°C in aerated water tanks (dissolved oxygen [DO] concentration: 7.8±0.5 mg/L) with a flow-through seawater supply. After 7 days of acclimation, hypoxia-exposure experiments were conducted at 25°C using published methods [3] by bubbling nitrogen gas into an aquarium. The desired concentration of DO was detected by using a DO meter (YSI, Canada). L. crocea can not maintain the aerobic pathway at DO levels below 2.0 mg/L, and it resorts to anaerobic metabolism [15]. Therefore, at the onset of hypoxia, the oxygen content in the tank was lowered from 7.8±0.5 mg/L to 1.6±0.2 mg/L over a 10-min period. Brains were harvested from six fish at the 0, 1, 3, 6, 12, 24, and 48 h time points and frozen immediately in liquid nitrogen until RNA extraction and transcriptome analyses were performed. To investigate whether the gene expression alterations under hypoxia come from changes in the baseline expression levels, brain tissues were collected from five fish with no hypoxia treatment at several corresponding time points (0, 3, 6, 12, and 24 h) and subjected to real-time PCR analysis of selected genes (S19 Fig).

Total RNA was extracted from the brain tissues harvested at 0, 1, 3, 6, 12, 24, and 48 h after hypoxia exposure using the guanidinium thiocyanate-phenol-chloroform extraction method (Trizol, Invitrogen, USA). The RNA samples (5 μg each) were used to construct RNA-seq libraries using the Illumina mRNA-Seq Prep Kit. The libraries were sequenced by using the Illumina HiSeq 2000 sequencing platform with the paired-end sequencing module [76]. More than 4 gigabase (Gb) of clean data were generated in each library (S29 Table). Clean reads were picked out from the raw reads of each library following removal of reads containing adaptor sequences, reads with an N (unknown bases in read) percentage higher than 5% and low quality reads (>30% of bases with a quantity score Q-value ≤ 10) using an in-house C++ program.

Clean reads were aligned to the L. crocea genome and its coding sequences (CDS) which were produced from genome-wide prediction of protein-coding genes (S29 Table), with SOAPaligner/SOAP2 [60]. The alignment to the CDS of L. crocea genome was then utilised to calculate the distribution of reads on reference genes and to perform coverage analysis. If an alignment result passed quality control (alignment ratio > 70%), we proceeded to gene expression calculations and differential expression comparisons. The gene expression levels were calculated based on RPKM values [69], and comparisons of gene expression difference between control (0 h) and each time point after hypoxia exposure (1, 3, 6, 12, 24, and 48 h) were performed using the method described by Audic and Claverie [77]. Genes with fold change ≥ 2 and FDR ≤ 0.001 were defined as differentially expressed genes (DEGs) [78]. Gene expression patterns were analyzed by the open source clustering software Cluster 3.0 and presented using the GraphPad Prism 5 software or Java TreeView. Raw sequencing data for the transcriptome have been deposited in the Gene Expression Omnibus (GEO) under accession GSE57608.

Real-time PCR analyses of differentially expressed genes under hypoxia

Real-time PCR analysis was performed using the Mastercycler epgradient realplex4 (Eppendorf, Germany) with SYBR Green as the fluorescent dye, according to the manufacturer’s protocol (Takara, China). Primer set was designed based on the coding sequence of each identified gene in L. crocea genome (S30 Table). Total RNA was extracted from brain tissues of five fish sampled at 0, 1, 3, 6, 12, 24, and 48 h after hypoxia induction. First-strand cDNA was synthesized from 2 μg total RNA and used as a template for real-time PCR. Real-time PCR was performed in a total volume of 20 μL, and cycling conditions were 95°C for 1 min, followed by 40 cycles of 94°C for 5 s, 57°C for 15 s, and 72°C for 20 s. The expression levels of each gene were expressed relative to those of β-actin in each sample using the 2^-ΔΔCT method [79]. Each real-time PCR assay was repeated three times. The data of real-time PCR were expressed as the standard error of the mean (SEM). Two-tailed Student’s t test was used for the significance test of the gene expression levels in brain tissues between 0 h and each time point after hypoxia exposure. A P-value <0.05 was considered to be statistically significant.

LC—MS/MS analyses and mucus protein identification

Skin mucus was collected from six healthy L. crocea individuals under air exposure as previously described [50]. Briefly, the fish were anesthetised with a sub-lethal dose of Tricaine-S (100 mg/L), and transferred gently to a sterile plastic bag for 3 min to slough off the mucus under air exposure. To exclude the cell contamination, mucus was diluted in fresh, cold phosphate-buffered saline and drop-splashed onto slides, which were then air-dried. After staining with 10% Giemsa dye (Sigma, St Louis, MO, USA) for 20 min, the mucus was observed under a Nikon microscope with a 20 × objective. No cell was observed.

Proteins were extracted from a pool of skin mucus of six fish by the trichloroacetic acid-acetone precipitation method and digested by the trypsin gold (Promega, USA). The peptides were then separated by the strong cation exchange chromatography using a Shimadzu LC-20AB HPLC Pump system (Kyoto, Japan). Data acquisition was performed with a Triple TOF 5600 System (AB SCIEX, Concord, ON) fitted with a Nanospray III source (AB SCIEX, Concord, ON). All spectra were mapped by MASCOT server version 2.3.02 against the database of the L. crocea annotated proteome with the parameters as follows: peptide mass tolerance 0.05 Da; fragment mass tolerance 0.1 Da; fixed modification “Carbamidomethyl (C)”; variable modifications “Gln->pyro-Glu (N-term Q), Oxidation (M), and Deamidation (N, Q)”. The ion score of a peptide matched to a protein must be greater than or equal to the MASCOT identity score, and the peptides must have a length of at least 6 amino acids. For further analysis of the function of the mucus proteome, only the proteins with at least two unique peptides were selected. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) with the dataset identifier PXD001218.