Draft genome of the protandrous Chinese black porgy, Acanthopagrus schlegelii.

Background: As one of the most popular and valuable commercial marine fishes in China and East Asian countries, the Chinese black porgy (Acanthopagrus schlegelii), also known as the blackhead seabream, has some attractive characteristics such as fast growth rate, good meat quality, resistance to diseases, and excellent adaptability to various environments. Furthermore, the black porgy is a good model for investigating sex changes in fish due to its protandrous hermaphroditism. Here, we obtained a high-quality genome assembly of this interesting teleost species and performed a genomic survey on potential genes associated with the sex-change phenomenon. Findings: We generated 175.4 gigabases (Gb) of clean sequence reads using a whole-genome shotgun sequencing strategy. The final genome assembly is approximately 688.1 megabases (Mb), accounting for 93% of the estimated genome size (739.6 Mb). The achieved scaffold N50 is 7.6 Mb, reaching a relatively high level among sequenced fish species. We identified 19 465 protein-coding genes, which had an average transcript length of 17.3 kb. By performing a comparative genomic analysis, we found 3 types of genes potentially associated with sex change, which are useful for studying the genetic basis of the protandrous hermaphroditism. Conclusions: We provide a draft genome assembly of the Chinese black porgy and discuss the potential genetic mechanisms of sex change. These data are also an important resource for studying the biology and for facilitating breeding of this economically important fish.

Data Description

Background information

As one of the most popular and valuable commercial marine fishes in China and East Asian countries, the Chinese black porgy (Acanthopagrus schlegelii), also known as the blackhead seabream, has some interesting characteristics such as fast growth rate, good meat quality, resistance to diseases, and good adaptability to various environments. It is often farmed for food in the South China Sea and the coastal waters around Japan and Korea [1,2]. In addition, it is an eurythermal and euryhaline fish, living in a wide range of water temperatures and salinities. Recently, some basic studies on the genetic improvement of its growth and its disease resistance have been performed in order to increase efficiency of farming [3].

The Chinese black porgy is also a good model for investigating the genetic mechanisms of sex change due to its interesting life cycle. It is a functional male during the first 2 years and subsequently a female during the next couple of years. Recently, a good hybrid of the Japanese seabream (Pagrosomus major; ♀) and the Chinese black porgy (♂) has become available [4,5], with better growth performance and higher tolerance against low temperature than its parents.

However, the genetic mechanisms for these interesting biological characteristics are still unclear. Here, we sequenced and assembled the whole genome of the Chinese black porgy, before performing a genomic survey on potential genes associated with the sex-change phenomenon.

Sample and sequencing

The wild black porgy (National Center for Biotechnology Information [NCBI] Taxonomy ID: 72 011; Fishbase ID: 6531) individuals (Fig. 1) were collected from Laizhou Bay in Yantai, Shandong Province, China. Genomic DNA was extracted from the muscle of a female fish using Qiagen GenomicTip100 (Qiagen, Hilden, USA). We used the whole-genome shotgun sequencing strategy and constructed the subsequent 3 short-insert libraries (250-bp, 500-bp, and 800-bp) and 4 long-insert libraries (2-kb, 5-kb, 10-kb, and 20-kb) in accordance with the standard protocol from Illumina (San Diego, USA). These constructed libraries were sequenced on the Illumina HiSeq 2000 system [6] (the read length is 125 bp). Finally, we generated 272.9-Gb raw reads from all 7 libraries.

Figure 1:

Image of a Chinese black porgy captured from Laizhou Bay in Yantai, Shandong Province, China.

Before assembly of the sequencing reads, SOAPfilter v2.2 software [7] with default parameters (-y -p -g 1 -o clean -M 2 -f 0) was used to remove low-quality raw reads (including reads with 10 or more nonsequenced/low-quality bases), polymerase chain reaction duplicates, and adaptor sequences. Subsequently, we obtained approximately 175.4 Gb of clean reads for further genome size prediction and assembling. A k-mer analysis with the formula G = k_num/k_depth [8] was performed to estimate the genome size of Chinese black porgy. In our current study, the achieved total number of k-mers and k_depth was 2.81 × 1010 and 38, respectively. Therefore, the genome size of Chinese black porgy is estimated to be 739.6 Mb. Based on this result, the retained reads were calculated to cover approximately 238-fold of the whole genome.

Assembly and evaluation

To obtain a genome assembly, we used SOAPdenovo2 v2.04.4 (SOAPdenovo2, RRID:SCR_014986) [9] with optimized parameters (pre-graph -K 27 -p 16 -d 1; contig –M 3; scaff -F -b 1.5 -p 16) to deal with these clean reads. In brief, the reads from short-insert libraries were applied for the contig assembly, before alignment of all the filtered reads onto the contigs for linking these contigs to generate scaffolds. GapCloser v1.12 (GapCloser, RRID:SCR_015026) [7] with default parameters was subsequently used to fill some intra-scaffold gaps in the local assembly, in which the reads were equipped with one end uniquely mapped to a contig and the other end located within a gap. Meanwhile, SSPACE (version 2.0) [10] with default parameters was used to obtain super scaffolds with the reads from the long-insert libraries (2-kb, 5-kb, 10-kb, and 20-kb). The final genome assembly was approximately 688.1 Mb, which accounts for 93.0% of the estimated genome size (739.6 Mb; Table 1).

Table 1:

Summary of the genome assembly and annotation

Genome assembly parameter
contig N50 size, kb	17.2
contig number, > 100 bp	115 091
Scaffold N50 size, Mb	7.6
Scaffold number, > 100 bp	31 359
Total length, Mb	688.1
Genome coverage, ×	257.6
Longest scaffold, bp	22 574 836
Genome annotation parameter
Protein-coding gene number	19 465
Mean transcript length, kb	17.3
Mean exons per gene	11.1
Mean exon length, bp	170.2
Mean intron length, bp	1519.2

The achieved scaffold N50 is 7.64 Mb, reaching a relatively high length among sequenced fish species. In comparison, other scaffolds have levels of 1.55 Mb for the zebrafish [11], 1.1 Mb for platy fish [12], 867 kb for half-smooth tongue sole [13], 1 Mb for common carp [14], 6.4 Mb for grass carp [15], 2.97 Mb for Atlantic salmon [16], 1.8 Mb for a seahorse [17], and 1.15 Mb for a Chinese barbel fish [18]. The Core Eukaryotic Genes Mapping Approach (CEGMA, RRID:SCR_015055), version 2.5 [19], with a set of 248 conserved core eukaryotic genes (CEGs) was used to assess the completeness of the final assembly. The estimates suggest that 90.7% CEGs are complete and 92.3% are partial. Meanwhile, Benchmarking Universal Single-Copy Orthologs (BUSCO, RRID:SCR_015008), version 3, [20] was applied to evaluate the quality of the generated genome assembly. We chose the representative actinopterygian gene set with 4584 single-copy genes as the reference. The BUSCO values were calculated as follows: C: 89.1% [S: 86.2%, D: 2.9%], F: 2.5%, M: 8.4%, n: 4584, in which percentages of the total gene number (n) for the complete (C), single (S), duplicated (D), fragmented (F), and missed (M) are clarified. The results from CEGMA and BUSCO suggest that the assembled genome covers the majority of the gene space.

Annotation

We used RepeatProteinMask (version 4.0.6) [21] in RepeatMasker (RepeatMasker, RRID:SCR_012954) to identify the repetitive sequences, before using RepeatModeller (version 1.05) [22] and LTR_FINDER.x86_64-1.0.6 to construct a de novo repeat library. Additionally, repetitive elements were predicted using Tandem Repeat Finder (version 4.04). Finally, we observed that the identified repeat sequences cover 19.78% of the assembled genome (Table 2).

Table 2:

Detailed classification of repeat sequences in the assembled genome

	Repbase TEs		TE proteins		De novo		Combined TEs
Type	Length, Mb	In genome, %	Length, Mb	In genome, %	Length, Mb	In genome, %	Length, Mb	In genome, %
DNA	20.930	3.041	2.200	0.320	58.340	8.479	68.130	9.902
LINE	10.240	1.488	6.950	1.010	26.760	3.889	33.020	4.789
SINE	1.120	0.163	2.340	0.000	3.780	0.550	4.550	0.661
LTR	7.200	1.046	35.410	0.340	25.980	3.062	31.270	4.544
Other	0.020	0.003	0.000	0.000	0.000	0.000	0.020	0.003
Unknown	0.000	0.000	0.000	0.000	25.370	3.687	25.370	3.687
Total	35.300	5.130	11.480	1.669	124.540	18.099	136.240	19.780

Prediction of protein-coding genes was performed based on the integration of ab initio prediction, homologue prediction, and transcriptome-based prediction. The ab initio prediction was carried out with Augustus (Augustus: Gene Prediction, RRID:SCR_008417), version 2.5 [23], and GENSCAN (GENSCAN, RRID:SCR_012902), version 1.0, [24], on the repeat-masked assembly. For the homology-based gene prediction, homologous proteins of several reported fishes (zebrafish, Japanese puffer, stickleback, and medaka) were downloaded from Ensembl release 75 and aligned to the assembled genome using tBlastn (version 2.2.19) with e-value ≤ 1e–5. Subsequently, all the achieved alignments were analyzed using Genewise (version 2.2.0) software [25] to search for precise gene structures. We further filtered out these short (less than 150 bp), prematurely terminated or frame-shifted genes. For the transcriptome-based prediction, we obtained transcriptome data from a mixture of liver, muscle, skin, gill, and brain of a female fish at cDNA level. Those with low-quality bases, adapter sequences, and duplicated sequences were removed, and we acquired approximately 8 Gb of high-quality clean reads. Subsequently, TopHat2.1.1 [26] and Cufflinks (Cufflinks, RRID:SCR_014597), version 2.2.1 [27], were applied to predict gene structures using these retained reads. Eventually, the 3 gene sets generated from the prediction approaches were integrated into a comprehensive and nonredundant gene set using GLEAN [28]. As summarized in Table 1, the final gene set contains 19 465 genes, with an average transcript length of 17.3 kb. In addition, we ran BUSCO v3 [20] on the predicted coding sequences (CDS), and the final BUSCO score was up to 85.5% (C:85.5% [S:82.3%, D:3.2%], F:2.8%, M:11.7%, n:4584).

Simultaneously, all the protein sequences from the GLEAN analysis were mapped onto several public databases, including Pfam [29], PRINTS [30], ProDom [31], and SMART [32], to detect the known motifs and domains within our genome assembly. The data demonstrated that 99.3% of the predicted genes from the assembled genome contain at least 1 related functional assignment from other public databases, including Swiss-Prot [33], Interpro [34], TrEMBL [35], and KEGG [36].

Phylogenetic analysis

In order to examine the phylogenetic position of the Chinese black porgy, we downloaded protein sequences of 7 reported fishes, including spotted gar (Lepisosteus oculatus), stickleback (Gasterosteus aculeatus), Japanses fugu (Takifugu rubripes), medaka (Oryzias latipes), zebrafish (Danio rerio), platyfish (Xiphophorus maculatus), and Nile tilapia (Oreochromis niloticus) from Ensembl (release 83) [37]. These sequences were used to construct gene families by OrthoMCL (OrthoMCL DB: Ortholog Groups of Protein Sequences, RRID:SCR_007839) [38] and eventually generated 17 431 gene families by the all-to-all Basic Local Alignment Search Tool for Proteins strategy with an e-value of 1e–5. In additional, 65 gene families were only presented in the black porgy genome.

Subsequently, 3239 single-copy orthologous genes from these gene families were selected. These single-copy genes were further aligned using MUSCLE (MUSCLE, RRID:SCR_011812), version 3.8.31, with default parameters [39], before the protein alignments were changed to corresponding CDS using an in-house perl script. All nucleotide sequences of each species were integrated into a supergene, which was used to build a phylogenetic tree with PhyML (PhyML, RRID:SCR_014629) [40]. Our final data orientated the phylogenetic position of the black porgy in teleost (Fig. 2).

Figure 2:

Phylogeny of ray-finned fishes. Spotted gar was used as the outgroup. The bootstrap support value for the topology is 100. The pictures in the phylogenetic tree were downloaded from Fishbase.

Analysis of 3 types of genes for sex change

Sex change (secondary sex determination) is a universal phenomenon in fish but it usually does not occur in amphibians and mammals. The black porgy is a good model for the study on the molecular mechanisms of sex change. For providing a genomic survey on these genes in the assembled genome, protein sequences of 3 main types of genes potentially associated with sex change, including sex determination and differentiation genes, pluripotency factors, and apoptosis factors [41-43], were downloaded from the NCBI database and used for homology searches against the black porgy genome with tBlastn (version 2.2.19) [44]. We chose alignments with coverage >70% and identity >70% for further prediction of gene structures using Genewise (version 2.2.0) [25]. Finally, we obtained homologous sequences of 26 genes in the genome assembly of Chinese black porgy (for more details, see Table 3). All these predicted protein sequences were used to build a phylogenetic tree using PhyML [40], and we eventually observed that they were clustered with each corresponding homologue from other vertebrates.

Table 3:

Three types of genes potentially related to sex change in the black porgy genome

Sex determination and differentiation genes
Gene	Copy number	Scaffold	Copy number in other teleosts
fst	2	10, 17	1 (zebrafish, medaka)
sox9	5	11, 13, 16, 19, 27	2 (zebrafish, catfish)
vasa	10	11, 14, 16, 20, 27, 34, 37, 47, 53, 68	1 (zebrafish), 3 (Nile tilapia)
ctnnb1	4	2, 16, 64, 115	1 (zebrafish)
piwil1	1	15	1 (zebrafish, medaka)
piwil2	1	15	1 (zebrafish, medaka)
sf1	1	108	1 (zebrafish, medaka)
rspo1	2	2, 74	1 (zebrafish, medaka)
foxl2	2	1, 22	1 (zebrafish, medaka)
cyp19a1a	2	8, 28	1 (zebrafish, medaka)
gsdf	1	3	1 (medaka)
figla	1	32	1 (zebrafish, medaka)
dmrt1	1	56	1 (zebrafish, medaka, tongue sole)
wnt4	15	1, 2, 5, 6, 7, 8, 9, 18, 19, 20, 32, 34, 62, 67, 122	2 (zebrafish), 3 (rainbow trout)
dax1	4	2, 3, 14, 43	1 (medaka, tongue sole)
cyp11a1	2	8, 33	1 (zebrafish)
hsd3b1	2	7, 36	1 (zebrafish, medaka)
amhr2	2	9, 185	1 (medaka)
jnk1	9	1, 3, 4, 5, 16, 17, 38, 79, 117	1 (zebrafish)
Pluripotency factors
klf4	5	1, 3, 17, 96, 142	1 (zebrafish, medaka)
nr5a2	3	8, 19, 28	1 (zebrafish, medaka)
lin28a	2	2, 3	1 (zebrafish)
pou2	1	3	1 (zebrafish)
Apoptosis factors
traf2	2	3, 15	1 (zebrafish, medaka)
casp2	1	2	1 (zebrafish)
tnfr1	1	2	1 (zebrafish, medaka, tilapia)

Note: The last column states the gene copy number in other teleosts based on the phylogenetic trees (uploaded to GigaDB [69]) in this study.

Previous studies have revealed that multiple genes, including dmrt1, cyp19a1a, wnt4, sox9, sf1, foxl2, figla, amhr2, and dax1, are associated with sex change in the black porgy [41,45-47]. These sex determination and differentiation genes were also identified in our assembled scaffolds (in the first batch in Table 3). In the current study, the important male-related dmrt1 and the steroidogenesis-suppressing factor dax1 were mapped on scaffolds 56 and 14 of the black porgy genome, respectively.

It was reported that dmrt1 may play a key role in the sex change of the black porgy, while the male-phase maintenance of male development was regulated by the brain–pituitary–gonadal axis via the GnRH-GtH-Dmrt1 pathway [41]. In the economically important half-smooth tongue sole (Cynoglossus semilaevis), dmrt1 has been proven to be a necessary male sex-determining gene [48,49]. Moreover, previous findings suggest that a duplicate of dmrt1 is the male sex determinant in medaka and dmrt1 mutation causes a male–female sex reversal [50,51]. We also validated the existence of foxl2 and cyp19a1a, 2 putative female-related genes, in the black porgy genome. Previous findings revealed that cyp19a1a plays dual roles in gonadal development, while both cyp19a1a and foxl2 are related to the sex change of the black porgy [47]. However, foxl2 has proved to participate in sex differentiation, although it is not essential for the sex determination and sex change in the tongue sole [52].

With only one copy in the black porgy, figla is a germ-cell–specific transcription factor related to ovary development and differentiation [53]. However, 2 isotypes (figla_tv1 and figla_tv2) were reported in the tongue sole. It is noted that figla_tv1 possesses a conserved function in folliculogenesis as found in other vertebrates, while figla_tv2 may play a role in the spermatogenesis of pseudo-males by regulating the synthesis and metabolism of steroid hormones [53]. Also identified with 1 gene in the black porgy, sf1 was reported to act as an essential transcriptional factor for steroidogenesis and for development of the reproductive axis (Table 3) [54].

Interestingly, 5 copies of sox9 were also identified in the black porgy genome. Nevertheless, previous findings reported that only 2 paralogs of sox9 (sox9a and sox9b) are present in zebrafish [55] and catfish [56]. Paralog sox9a is usually associated with testicular development, while this may be linked with sex reversal in the tong sole [52]. In comparison, sox9b possesses a new function in the ovary [55]. In addition, we noticed that female-related genes (wnt4, vasa, and jnk1) had multiple copies in our current study, which may be retained since the whole-genome duplication in the ancestor of the teleost (Table 3). These genes have been proven to play important roles in ovarian growth and natural sex changes in fishes [57-60]. It was reported that 2 wnt4 genes (wnt4a and wnt4b) are present in most teleost fish, while other vertebrates and invertebrates possess only a single wnt4 gene. Furthermore, 2 copies of wnt4a, wnt4a1, and wnt4a2 exist in some teleost species resulting from the additional duplication of the wnt4 gene [61]. It has been shown that wnt4a was mainly expressed in the gonad, gill, and brain of teleost fish (such as zebrafish [62] and rainbow trout [63]), and it was confirmed to be associated with sex reversal in the tongue sole [61]. The vasa gene, also called ddx4, was reported to play an important role in gametogenesis and germ cell development [64]. Previous findings showed that vasa was a single copy gene in the majority of chordates such as zebrafish [65,66]. However, 3 vasa genes were also reported in Nile tilapia (Oreochromis niloticus) [67]. Jnk1 is closely associated with ovarian differentiation and development in fish. A previous finding [58] reported that jnk1 highly transcribed in the ovary of the female ricefield eel (Monopterus albus), another teleost with natural sex-change from female to male, and reduced to a substantial level at the subsequent stage of intersex. Hence, the data demonstrated that jnk1 may play a key role in sexual reversal. Surprisingly, 2 jnk1 genes (jnk1a and jnk1b) were reported in the polyploid hybrids of red crucian carp (Carassius auratus red var.) and common carp (Cyprinus carpio L.) [68].

Interestingly, our data demonstrate that the distribution of these 3 types of genes in the black porgy genome is similar to that in ricefield eel (our unpublished results; data from the Monopterus Whole Genome Shotgun Project have been deposited at DDBJ/EMBL/GenBank under accession number AONE00000000). For example, 2 male-related genes (piwil1 and piwil2) are clustered together, while lin28a and rspo1 are adjacent to each other. We also observed that most of these genes are congregated on scaffolds 1, 2, 3, 11, and 15 (Table 3).

Conclusions

In summary, we sequenced and assembled the whole genome of Chinese black porgy. This is the first genomic report of Sparidae fish. Furthermore, we provided a genomic survey on the 26 genes potentially associated with sex change. The achieved genome data will be helpful for further biological and evolutionary studies. Furthermore, it will be valuable for implementation of molecular breeding, with substantial support from our genomic data, to obtain genetic improvement of this economically important teleost fish.

Ethics approval and consent to participate

All animal experiments in this study were implemented according to the guidelines of the Animal Ethics Committee and ratified by the Institutional Review Board of Bioethics and Biosafety of BGI, China.

Availability of supporting data

The raw sequencing reads of all libraries and the transcriptome data have been deposited in the NCBI SRA database (accession numbers SRA541936 and SRA587358). Supporting data are available in the GigaScience database, GigaDB [69].

Competing interests

All authors report no competing interests.

Abbreviations

BUSCO: Benchmarking Universal Single-Copy Orthologs; CDS: coding sequence; CEGMA: Core Eukaryotic Genes Mapping Approach; CEG: core eukaryotic gene; NCBI, National Center for Biotechnology Information

Author contributions

Zy.Z., Q.S., and P.X. conceived the project. J.X., C.J., J.Q., F.Z., Hx.L., Hl.L., D.S., Z.R., and J.C. extracted the genomic DNA and performed genome sequencing. K.Z., S.C., Zw.Z., X.Y., J.Z., C.B., and J.L. assembled the genome and analyzed the data. T.G., R.G., and J.X. participated in discussions and provided valuable advice for revision. K.Z., Q.S., Zy.Z., P.X., Zw.Z., and S.C. prepared the manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

05 Sep 2017 Reviewed

Click here for additional data file.

07 Nov 2017 Reviewed

Click here for additional data file.

11 Sep 2017 Reviewed

Click here for additional data file.

07 Nov 2017 Reviewed

Click here for additional data file.

08 Dec 2017 Reviewed

Click here for additional data file.