Genome-wide SNP identification for the construction of a high-resolution genetic map of Japanese flounder (Paralichthys olivaceus): applications to QTL mapping of Vibrio anguillarum disease resistance and comparative genomic analysis.

High-resolution genetic maps are essential for fine mapping of complex traits, genome assembly, and comparative genomic analysis. Single-nucleotide polymorphisms (SNPs) are the primary molecular markers used for genetic map construction. In this study, we identified 13,362 SNPs evenly distributed across the Japanese flounder (Paralichthys olivaceus) genome. Of these SNPs, 12,712 high-confidence SNPs were subjected to high-throughput genotyping and assigned to 24 consensus linkage groups (LGs). The total length of the genetic linkage map was 3,497.29 cM with an average distance of 0.47 cM between loci, thereby representing the densest genetic map currently reported for Japanese flounder. Nine positive quantitative trait loci (QTLs) forming two main clusters for Vibrio anguillarum disease resistance were detected. All QTLs could explain 5.1-8.38% of the total phenotypic variation. Synteny analysis of the QTL regions on the genome assembly revealed 12 immune-related genes, among them 4 genes strongly associated with V. anguillarum disease resistance. In addition, 246 genome assembly scaffolds with an average size of 21.79 Mb were anchored onto the LGs; these scaffolds, comprising 522.99 Mb, represented 95.78% of assembled genomic sequences. The mapped assembly scaffolds in Japanese flounder were used for genome synteny analyses against zebrafish (Danio rerio) and medaka (Oryzias latipes). Flounder and medaka were found to possess almost one-to-one synteny, whereas flounder and zebrafish exhibited a multi-syntenic correspondence. The newly developed high-resolution genetic map, which will facilitate QTL mapping, scaffold assembly, and genome synteny analysis of Japanese flounder, marks a milestone in the ongoing genome project for this species.

Introduction

Japanese flounder (Paralichthys olivaceus), one of the most desirable and highly priced marine fish species, is widely cultured along the coast of Northeast Asian countries such as China, Japan, and Korea.[1] In addition to a booming Chinese market, >60,000 tons of Japanese flounder were artificially produced in Japan and Korea in 2013,[2] making P. olivaceus one of the most important species in stock enhancement programs in these countries. As a result of long-term, high-intensity stocking and resource management, however, farmed Japanese flounder have depressed immune systems that enhance individual susceptibility to microbial infections.[3] To increase profitability and sustainability while maintaining genetic variability in the cultured stock, the development of genetic breeding programs, such as marker- and gene-assisted selection, is urgently required. These approaches involve selection of genomic loci or genes related to the economic trait of interest using genetic maps constructed by molecular markers.[4]

An accurate genetic linkage map is therefore foundational to genetic breeding of a species.[5] To date, linkage maps have been constructed for at least 28 fish species; in addition, economically important traits have been mapped using these maps in over 13 fish species, including tilapia, rainbow trout, channel catfish, Atlantic salmon, common carp, and Asian seabass.[6] Following the construction of the first generation of Japanese flounder genetic maps based on simple sequence repeat (SSR) and amplified fragment length polymorphism (AFLP) markers by Coimbra et al.[7] and SSRs alone by Kang et al.,[8] the second generation of linkage maps has been developed with various types of molecular markers, such as SSRs and a few single-nucleotide polymorphism (SNP) markers, either alone or in combination.[9,10] Such efforts have resulted in the identification of several quantitative trait loci (QTLs) for disease resistance and growth in Japanese flounder. For instance, a single major genetic locus associated with lymphocystis disease resistance has been detected[11] and its candidate gene tlr2 identified,[12] making the development of a lymphocystis disease-resistant strain of Japanese flounder possible by marker-assisted selection.[13] In addition, four QTLs for growth rate are clustered on genetic linkage group (LG) 14, suggesting that these major loci contributing to growth trait variation are potentially useful for marker-assisted selection in future genetic breeding programs of Japanese flounder.[10] Although previous genetic maps of Japanese flounder have been successfully used to map QTLs for a few economic traits, the number of commonly available markers, restricted to hundreds to a few thousand, makes it difficult to carry out fine-scale mapping to shrink the genomic regions tightly associated with important traits. A genetic map of much higher resolution is thus urgently needed for Japanese flounder.

A high-resolution genetic linkage map is also an excellent tool for genome map construction, which in turn allows direct comparison of chromosomal organization and evolution.[14] For instance, Caenorhabditis briggsae chromosomes reconstructed by SNP-based genetic mapping were compared with the Caenorhabditis elegans genome. An almost complete conservation of synteny was observed, with many cases of perfect 1:1 orthologues found between these two distantly related species.[14] In addition, construction of a dense gene-based map of the butterfly species Bicyclus anynana enabled broad-coverage analysis of synteny with the lepidopteran reference genome; this analysis suggested strong conservation of gene assignments to chromosomes and numerous large- and small-scale chromosomal rearrangements.[15] The application of genetic mapping allows genome-level analysis to be extended to non-model species, conferring an advantage in the comparative analysis of chromosomal organization and evolution.

With the emergence of next-generation (massively parallel) sequencing and associated technological genotyping advancements such as reduced representation library sequencing techniques,[16] SNP markers, representing the most abundant source of variation in the genome, have become widely used for high-resolution genetic map construction; they are also applied for measurement of genetic diversity and investigation of population structure.[17] As a reliable, high-throughput, affordable method to reduce genomic complexity, restriction-site associated DNA tag sequencing (RAD-seq) has been particularly attractive for SNP discovery and genotyping.[18] RAD-seq technology has been successfully applied in the construction of genetic maps for various species, such as stickleback,[19] rainbow trout,[20] eggplant,[21] barley,[22] chickpea,[23] guppy,[24] and lupin.[25] Furthermore, the discovery that replicate parallel phenotypic evolution may be occurring in stickleback on a genome-wide scale[26] as well as the detection of population admixture and improved identification of hybrid and non-hybridized individuals using RAD-seq[20] suggest that this technique is applicable to population genomics and phylogeography.

Despite the availability of two generations of genetic maps, the large-scale discovery and utilization of SNPs, the most promising markers for providing sufficiently dense genetic maps,[23] has not been carried out in Japanese flounder. In the present study, we consequently performed a large-scale identification of genome-wide SNPs derived from RAD-seq of a Japanese flounder mapping population containing 2 parents and 216 offspring. We then used the results to construct a high-density (third generation) SNP-based map. We also successfully detected QTLs for Vibrio anguillarum disease resistance and their related genes in Japanese flounder. Finally, we used our SNP-based genetic map in conjunction with the P. olivaceus Genome Sequencing Project (PoGSP) database, which is to be published in the near future. In particular, the map was used to facilitate the anchoring and orienting of scaffolds generated by whole genome sequence data assembly and for chromosomal-level comparative analysis.

Materials and methods

Mapping population

Twenty-four full-sib families of Japanese flounder were bred by Haiyang Yellow Sea Fisheries Co. (Yantai, China) in 2010.[27] All 216 offspring of one of the families, with an average length of 12.61 ± 1.40 cm and average weight of 18.85 ± 6.28 g, were challenged through intraperitoneal injection with a 0.2-ml bacterial suspension of ∼7.16 × 105 colony-forming units (CFU) of V. anguillarum. Another 20 individuals were injected with 1× phosphate-buffered saline as a control. The median lethal concentration had been determined in our laboratory previously.[28] The offspring were kept in three 0.28-m3 tanks supplied with fresh seawater at 18 ± 0.5°C. This challenge experiment was performed once and lasted for ∼1 week. Mortality was recorded every 3 h based on the appearance of dead fish (Supplementary Table S1). Genomic DNA was isolated from the fins of the parental fish and 216 offspring using traditional phenol–chloroform extraction in combination with RNase treatment and stored at −20°C.[10] Before construction of RAD-seq libraries, all DNA samples were quantified using a NanoDrop instrument (Thermo Scientific, Wilmington, DE, USA), and their concentrations were adjusted to 50 ng/µl using Tris–EDTA buffer.

RAD library construction and sequencing

RAD-seq libraries were constructed using a protocol adapted from Baird et al.[19] Briefly, 250 ng of genomic DNA from each of 218 library individuals was digested separately with 20 units of PstI and then heat inactivated at 65°C. Various P1 adapters, each with a unique 4–8 bp molecular-identifying sequence (MID), were then ligated to designated individuals, which were then pooled in groups of 24 individuals and randomly sheared to DNA fragments. Sheared DNA was purified, eluted, and separated using gel electrophoresis, and a DNA fraction corresponding to 300–700 bp was excised and purified. After end repair, purification, and elution, dATP overhangs were added to the DNA fraction. A paired-end P2 adapter containing T overhangs was ligated to 20 μl of sheared, size-selected, P1-ligated, and pooled DNA template with a specific adapter. The ligated material was then purified, eluted, and subjected to PCR enrichment. Sequencing of the RAD products from the 218 individuals was performed on a HiSeq2000 next-generation sequencing platform. Sequencing data for each individual were then extracted according to the specific MID.

SNP discovery and genotyping

We first filtered out Illumina short reads lacking sample-specific MIDs and expected restriction enzyme motifs. All the short reads from each of the samples were then clustered into tag reads on the basis of sequence similarity (allowing one mismatch at most between any two reads within each tag read cluster, with clusters having <2 or >100 reads discarded).[29] Tag reads from the two parents were compared and filtered to remove monomorphic DNA sequences, leaving only the tag reads with SNPs. SNPs were selected according to the following criteria: (i) alleles with a minimum coverage of five reads and a score >20 (P > 0.05) were selected; (ii) the base at the SNP site was unique (if at a homozygous locus) or bimorphic (if at a heterozygous locus, with the minor base represented by at least three reads); and (iii) the ratio of the two kinds of bases (major to minor) at a heterozygous locus ranged from 1 to 5. Regions containing these putative SNPs were used as reference SNP regions. The clean tag reads of offspring individuals were then aligned to reference SNP regions, with genotypes of individuals determined by reference to the parental genotypes. Filtering of reads and SNP selection were performed following previously described algorithms[29,30] as implemented using custom scripts (available from the authors upon request).

Genetic map construction

A pseudo-testcross population was used to construct the linkage map. For the linkage analysis, RAD-based SNPs were first tested against the expected segregation ratio. Two SNP alleles heterozygous with respect to the two parents were expected to segregate in a 1:2:1 ratio, whereas a pair of SNP alleles, one heterozygous and the other homozygous, were expected to segregate in a 1:1 ratio. To construct a genetic map, markers showing significant segregation distortion (P < 0.01, χ2 test) were removed. The remaining SNPs were then used to construct the genetic map with the recently developed software package Lep-MAP.[31] The genotype data were first filtered manually to remove obvious Mendelian errors from the offspring. LG assignments were then obtained using the separate chromosomes module with a logarithm of odds (LOD) score limit of 10. The marker order was obtained using the order markers module of Lep-MAP. Because this module uses only paternally informative markers in the ordering, a second genotype data file was constructed by swapping the parents of the original data. An integrated map was constructed by using both genotype files as input, with paternal and maternal maps generated by using only one input file. To speed up the computation, constant rates for genotype errors and recombinations were used. Finally, the marker positions and error parameters were established by order markers module.

QTL mapping for V. anguillarum disease resistance

QTL analyses were conducted with the WinQTLCart2.5 software program using the composite interval mapping (CIM) method.[32] The CIM analysis was run using Model 6 with four parameters for forward and backward stepwise regression, a 10-cM window size, five control markers, and a 1-cM step size. The LOD threshold value was determined on the basis of 1,000 permutations at a whole genome-wide significance level of P < 0.05. The location of each QTL was determined according to its LOD peak location and surrounding region. Other QTL reference values, including phenotypic variation and a positive or negative additive effect for V. anguillarum disease resistance, were also calculated by WinQTLCart2.5. Candidate genes for V. anguillarum disease resistance were identified by mapping the corresponding tags of SNPs in QTL regions to the scaffold assembly followed by retrieval of the corresponding gene ID from the gene annotation file.

Genome scaffold assembly and genome synteny

We used the high-resolution RAD-based SNP genetic map and BLASTN (E-value < 1 × 10−5, identity ≥95%, and alignment rate >50%) to map reads containing SNPs with less than two mismatches to PoGSP-derived scaffolds. To reduce the complexity of scaffolds anchored to hundreds of SNP markers, a tag SNP was selected from each scaffold with multiple SNPs. The position and orientation of scaffolds were ordered and assembled onto 24 pseudo-chromosomes corresponding to the 24 LGs based on genetic distances between SNPs. Scaffolds with only one SNP marker could be anchored but not oriented because of the lack of markers. We then used the Japanese flounder reference gene set generated from the PoGSP database along with medaka (Oryzias latipes) and zebrafish (Danio rerio) gene sequences from Ensembl (release 57). After filtering short genes (coding sequence <150 bp), we chose the transcripts with the longest coding sequences to represent each gene and then identified reciprocal best-matching orthologous genes between Japanese flounder and the other fish using BLASTP (E-value = 1 × 10−10). Syntenic regions defined by the top hits of the homology search were plotted using an in-house script.

Results

RAD-seq library construction and sequencing

A total of 218 RAD-seq libraries from 2 parents and their 216 offspring were constructed and sequenced on an Illumina HiSeq2000 platform to generate 3.91 billion raw reads. After data trimming, 3.52 billion reads, comprising ∼149.87 Gb of sequencing data, were individually partitioned into RAD tags according to their MIDs. Finally, female and male parental data sets, containing respectively 22.93 million filtered reads (comprising 947.96 Mb of data with a GC% of 47.02) and 17.7 million filtered reads (comprising 796.52 Mb of data with a GC% of 47.21), were correspondingly partitioned into 6,294,857 and 4,843,594 RAD tags. From the 216 offspring, a total of 3.48 billion filtered reads (average of 16.12 million) corresponding to 148,128.63 Mb of data (average of 685.78 Mb) were produced and divided into 939,770,718 RAD tags (ranging from 1,627,621 to 8,329,563 with an average of 4,350,790) for individual SNP discovery (Supplementary Table S1).

After stringent selection according to the above-described method, RAD tags from each individual were clustered and compared. RAD tags containing SNPs of the two parents were identified using a custom k-mer matching algorithm that excluded exact sequence matches (monomorphic loci) per 41-bp sequence. A panel of 13,362 high-fidelity SNPs with fixed genotypes in both parents was identified using these criteria described under ‘Material and methods’, and alleles for each marker were assigned to their respective parental donor. Using an in-house script, these SNPs were then analysed using genotyping data across the 216 offspring and classified into three categories: maternal heterozygous (6,741 SNPs), paternal heterozygous (4,128 SNPs), and heterozygous in both (2,493 SNPs). All 13,362 SNPs and their flanking sequences are listed in Supplementary Table S2.

High-resolution genetic map construction

A high-resolution RAD-based SNP genetic map of Japanese flounder based on a pseudo-testcross population was first constructed using Lep-MAP. A total of 12,712 segregating SNPs (95.14%) were successfully classified into 24 LGs (Table 1 and Fig. 1). The paternal map contained 8,784 SNPs with a total genetic distance of 2,561.17 cM; the length of each LG ranged from 19.53 cM (LG8) to 145.08 cM (LG19), with an average genetic length of 106.72 cM (Supplementary Fig. S1). The corresponding maternal map consisted of 6,343 SNPs representing a total length of 2,356.86 cM and ranging from 0.6 cM (LG8) to 187.33 cM (LG1) (Supplementary Fig. S2). More interestingly, LG8 was the smallest LG on both maternal (0.6 cM) and paternal maps (19.53 cM). LG8 on the maternal map contained only 3 SNPs heterozygous in both parents; on the paternal map, LG8 included 135 paternal heterozygous SNPs and 3 SNPs heterozygous in both parents. The resulting integrated map consisted of 24 LGs including 12,712 SNPs, which corresponded to 7,430 effective loci. The total map length was 3,497.29 cM, with an average inter-locus distance of 0.47 cM; the genetic length of each LG ranged from 21.64 cM (LG8) to 221.14 cM (LG20), with an average inter-locus distance of 0.37–0.75 cM. LG1 was the densest, having 495 effective loci with an average density of 0.42 cM, whereas LG8 had the least number of effective loci (only 59). On average, each LG contained 310 effective loci spanning 145.72 cM (Table 1 and Fig. 1). Locus names and SNP positions on the 24 LGs of the integrated genetic map are listed in Supplementary Table S3.

Table 1.

Characteristics of genetic maps and anchoring scaffolds of Japanese flounder

LG_ID	Paternal map		Maternal map		Integrated_Map
	No. of SNPs	Distance (cM)	No. of SNPs	Distance (cM)	No. of SNPs	No. of effective loci	Distance (cM)	Average inter-loci distance	No. of anchored scaffolds	Length of anchored scaffolds (Mb)
1	493	137.15	561	187.33	875	495	208.74	0.42	16	27.81
2	75	46.95	94	9.03	127	72	36.33	0.50	9	18.60
3	451	111.49	375	137.75	706	407	190.32	0.47	10	23.21
4	321	97.84	284	95.24	527	290	142.28	0.49	7	18.98
5	413	101.92	216	80.38	584	303	114.34	0.38	10	25.85
6	441	126.54	442	181.24	742	440	203.21	0.46	9	24.62
7	224	94.54	76	72.64	253	178	94.76	0.53	13	29.55
8	138	19.53	3	0.60	138	59	21.64	0.37	8	11.10
9	377	103.89	46	9.00	407	195	101.09	0.52	19	14.86
10	222	88.72	151	76.56	333	214	123.10	0.58	12	21.85
11	263	82.11	385	127.91	573	351	177.62	0.51	8	22.77
12	452	102.60	461	121.51	787	436	161.24	0.37	10	25.65
13	410	130.09	412	107.02	686	401	172.96	0.43	5	19.78
14	313	117.59	260	88.67	421	276	128.28	0.46	11	11.15
15	486	122.12	151	91.38	509	297	150.71	0.51	10	23.80
16	175	72.03	294	100.37	402	253	119.09	0.47	12	17.16
17	408	136.93	492	169.11	717	462	209.14	0.45	11	25.19
18	276	116.96	99	84.30	287	192	144.18	0.75	8	19.41
19	458	145.08	174	97.10	484	302	169.02	0.56	12	26.39
20	528	114.33	464	159.35	824	467	221.14	0.47	6	22.92
21	449	127.47	103	32.21	492	247	101.46	0.41	6	25.31
22	448	102.11	260	106.97	590	380	157.75	0.42	10	19.99
23	514	141.20	220	101.53	603	336	173.14	0.52	12	24.25
24	449	121.98	320	119.66	645	377	175.75	0.47	12	22.80
Total	8,784	2,561.17	6,343	2,356.86	12,712	7,430	3,497.29	0.47	246	522.99
Average	366	106.72	264	98.20	530	310	145.72	0.47	10	21.79

Figure 1.

Linkage group lengths and marker distribution of the high-resolution restriction site-associated DNA sequencing-based SNP genetic map of Japanese flounder. Within each linkage group, red, blue, and yellow lines, respectively, represent maternal heterozygous SNPs, paternal heterozygous SNPs, and SNPs heterozygous in both parents. Genetic map details are given in Supplementary Table S3.

Vibrio anguillarum disease resistance-associated QTLs and related genes

In total, 10 significant QTLs for V. anguillarum disease resistance were distributed on LG6, LG19, and LG21 of Japanese flounder (Table 2 and Fig. 2). Most of these QTLs were clustered together on their respective LGs. One major cluster containing six QTLs (qVA-1, qVA-2, qVA-3, qVA-4, qVA-5, and qVA-6) was detected between the narrow positions of 89.8–108.7 cM on LG6. Among them, qVA-2 located at 95.9–99.3 cM had the highest LOD value, 5, and correspondingly had the highest contribution to phenotypic variation, 8.38%. The other QTLs on LG6 were detected at positions 89.8–90.8, 100.7–102.5, 104.1–105.0, 105.0–105.4, and 107.8–108.7 cM, with LOD values of 3.0–3.6 and contributions to phenotypic variation of 5.1–7.59% (Supplementary Fig. S3). On LG21, another cluster situated within a short region (98.7–101.3 cM) consisted of three QTLs (qVA-8, qVA-9, and qVA-10) with a LOD value of 3.3–3.7 and was able to explain 6.06–6.56% of the phenotypic variation (Supplementary Fig. S3). Finally, qVA-7 containing a single SNP was centred around 115.6 cM on LG 19; it had a LOD value of 15.8 and explained 1.19% of the phenotypic variation. Although these 10 QTLs explained >59.25% of the total phenotypic variation, no major loci (explaining >20% of the total variation) were detected. The fact that these loci do not independently have higher contributions to such a complicated trait is not unexpected. Several genes related to V. anguillarum disease resistance were identified from the higher density map based on the Japanese flounder genome assembly. As shown in Supplementary Table S4, 218 genes were identified from the QTL regions, 12 of which were immune-related genes (tap1, rftn1, satb1, cd40, cd69, aicda, mtss1, ccr4, azi2, mrc1, nod1, and tgfbr2) functioning as key factors in different immune-gene signalling pathways or as antigen receptors. Interestingly, tap1 and satb1 act as molecular organizers for the major histocompatibility complex (MHC) class I,[33] while cd40 and cd69 are accessory molecules of MHC class II[34] directly involved in V. anguillarum disease resistance in Japanese flounder.

Table 2.

Characteristics of Vibrio anguillarum disease resistance QTLs

QTL	Linkage group	Genetic position	Associated marker	LOD	aExp%	Additive effect
qVA-1	LG6	89.8–90.8	record_231777.7	3	5.10	0.21
qVA-2	LG6	95.9–99.3	record_245752.29	5	8.38	0.31
qVA-3	LG6	100.7–102.5	record_247030.22	3.6	7.59	0.29
qVA-4	LG6	104.1–105	record_246301.16	3.3	5.87	0.26
qVA-5	LG6	105–105.4	record_254077.24	3.5	6.19	0.26
qVA-6	LG6	107.8–108.7	record_254548.37	3.4	5.95	0.26
qVA-7	LG19	115.6–115.8	record_255627.24	15.8	1.19	−0.12
qVA-8	LG21	98.7–99.6	record_245734.7	3.7	6.06	−0.25
qVA-9	LG21	100.6–100.9	record_255692.21	3.7	6.56	−0.27
qVA-10	LG21	100.9–101.3	record_245651.14	3.3	6.36	−0.26

aExp, percentage of explained phenotypic variation.

Figure 2.

Genetic location of QTLs for Vibrio anguillarum disease resistance along the Japanese flounder genome. The blue horizontal line represents a linkage group-wise logarithm of odds (LOD) significance threshold of 3.0.

Genome assembly and synteny

The high-resolution genetic map was used in the present study for genome scaffold assembly. Among the 12,712 SNPs on the high-resolution map, 12,463 were successfully used to anchor 246 scaffolds representing 522.99 Mb to create a genome map of Japanese flounder. In total, we constructed 24 pseudo-chromosomes based on the LGs; each chromosome comprised an average of 21.79 Mb (10 scaffolds), corresponding to a mean linkage distance of 145.72 cM and suggesting a physical/genetic distance ratio of 149.53 Kb/cM. Among the LGs, LG7 holding 13 scaffolds corresponded to the largest chromosome (29.55 Mb) and represented 95.78% of the entire assembly, whereas the smallest LGs contained 8 scaffolds representing 67.85% of the assembly (16.36 Mb) (Table 1). In addition, 219 scaffolds covering 520.44 Mb were oriented by assigning multiple genetically separated markers located on each scaffold. The orientation of 27 scaffolds covering 2.55 Mb could not be determined, because only single markers were present on those scaffolds.

Following genome map construction, we compared the chromosomal orders of protein-coding genes of Japanese flounder on each scaffold with their counterparts in medaka and zebrafish. A total of 6,768 1:1 best orthologues between Japanese flounder and zebrafish were identified. As shown in Supplementary Table S5 and Fig. 3, 7 of the 24 chromosomes of flounder (Po. 4, 6, 10, 11, 12, 15, 18, 21, and 22) were in relatively conserved synteny with zebrafish chromosomes Dr. 12, 16, 3, 2, 14, 9, 24, 15, and 19 respectively. Many of the flounder chromosomes were found to have syntenic blocks with more than one chromosome of zebrafish and vice versa. An example of the multi-syntenic correspondence between zebrafish and flounder genomes is that of Po. 20, which exhibited one of the highest degrees of synteny: Po. 20 had 519 hits on zebrafish chromosomes, most of which (171, 142, and 136) were located on Dr. 05, 10, and 21, respectively. Similarly, Dr. 06 was found to be highly syntenic with Po. 01, 05, and 24, whereas Po. 01, 05, and 24 were in synteny with Dr. 11, 01, and 22, respectively. With respect to flounder–medaka genome synteny, the alignment of 6,463 orthologues on the 24 flounder chromosomes revealed an obvious syntenic relationship with the 24 corresponding chromosomes of medaka. Except for some chromosomes with minor multi-chromosome hits (<50 genes), the results of this comparison suggest that flounder is more closely related to medaka than to zebrafish (Fig. 3).

Figure 3.

Circos diagram representing syntenic relationships between flounder and (a) zebrafish and (b) medaka, respectively. Each colored arc represents an orthologous match between two species. Pol, Paralichthys olivaceus; Dre, Danio rerio; Ola, Oryzias latipes.

Discussion

High-resolution genetic map construction using RAD-based SNPs

Genetic maps constructed using molecular markers are important for genomic and genetic analyses of individual species. SNPs are particularly attractive for genetic map construction, because they represent the most common type of DNA polymorphism in the genome and are amenable to high-throughput genotyping.[17] In the present study, ∼950 million RAD tags were generated from 2 Japanese flounder parents and 216 offspring individuals. We discovered 13,362 novel SNPs evenly distributed throughout the entire Japanese flounder genome, providing thereby a large number of genetic variation resources for future genome selection and genome-wide association studies. Among the identified SNPs, 12,712 SNPs were successfully genotyped and assigned to 24 LGs corresponding to the 24 chromosome pairs of the Japanese flounder genome. The total length of the resulting map was 3,497.29 cM, with an average marker spacing of 0.47 cM. Four genetic maps of Japanese flounder have been previously constructed using various markers,[7-10] but those maps are characterized by relatively larger genetic intervals and/or non-uniform marker distributions. The most recent genetic map[10] containing 1,487 SSRs with an average interval of 1.22 cM was formerly the densest reported flatfish linkage map and worthy of acclaim. However, whole genome-level analyses such as genome-wide association studies were still not practical using this map because of the limits imposed by the types of markers and their associated genotyping methods.[23] The denser SNP-based genetic map constructed in this study thus has greatly expanded application. It represents a considerable improvement over previous flounder genetic linkage maps based on SSR and AFLP markers.

Most LGs on our genetic map were considered to be saturated with an even distribution of three types of SNPs, namely, maternal heterozygous, paternal heterozygous, and double heterozygous SNPs; a few LGs, such as LG8, LG9, LG11, LG18, LG19, were not, however, revealing a complex genome containing regions enriched in repeat sequences and transposable elements.[35] Interestingly, LG8 encompassed the smallest genetic distance, contained only paternal and double heterozygous markers, and corresponded to the smallest physical distance (11.1 Mb). The fact that LG8 on the maternal map had only markers heterozygous in both parents whereas the group on the paternal map had only paternal heterozygous and double heterozygous markers suggests specific suppression of homologous recombination in this LG. Similar phenomena have been documented in other species having recently evolved sex chromosomes with large pseudoautosomal regions.[36,37] To overcome this complication, additional markers must be developed and a larger mapping population is required to increase the density of the linkage map.

The Lep-MAP program used in this study offers several advantages over JoinMap,[38] currently the most popular linkage-mapping software package. Lep-MAP, unlike JoinMAP, is suitable for use with tens of thousands of markers. Furthermore, Lep-MAP is faster on large data sets, requires no manual input, and handles genotyping errors better because it estimates genotyping error rates for each marker. Since its development in 2013, Lep-MAP has thus been heavily used for genetic map construction, such as for the high-resolution genetic map of Glanville fritillary.[39] Finally, the RAD-seq technique is proving to be a highly valuable tool for SNP discovery and genotyping followed by genetic map construction and has currently been applied to various species. The success of RAD-seq is owing to advances in next-generation technologies, which reduce the cost of DNA sequencing. Genotyping by sequencing is therefore expected to become increasingly popular for high-throughput genotyping, genetic map construction, and genome analysis.[40]

V. anguillarum disease resistance-related QTLs and associated genes

The generated high-resolution genetic map allowed us to perform QTL fine mapping for an economically important trait of Japanese flounder. In the present study, 10 QTLs associated with V. anguillarum disease resistance were found to be distributed on three LGs (LG6, LG19, and LG21). Interestingly, most of the QTLs were concentrated within a narrow region (cluster) on the LGs. Six QTLs were clustered together (89.8–108.7 cM) on LG6 corresponding to assembly scaffolds, scaffold196 (1.08 Mb), scaffold103 (1.06 Mb), scaffold31 (4.39 Mb), and scaffold36 (11.37 Mb), and three were found in another cluster (98.7–101.3 cM) on LG21, corresponding to scaffold190 (0.97 Mb) and scaffold70 (2.59 Mb). The quite small genetic and physical distances among QTLs in certain clusters suggest that the individual clusters are highly effective QTLs. One exception is qVA-7 on LG19. Even though it had the largest LOD value (15.8), qVA-7 contained only one SNP located in a physical gap between two scaffolds, indicating that it may be a false-positive QTL or, alternatively, the result of an error arising during genetic map construction due to segregation distortion. In addition, the smallest significant QTL (other than qVA-7) explained 5.1% of the phenotypic variation, whereas the maximum contribution was only 8.38%, reflecting the complexity of this polygenically controlled disease trait in Japanese flounder. Although the possibility exists that additional disease-controlling QTLs are present, the positive QTLs detected in this study, taken together, explained 58.06% of the total disease resistance variation, suggesting that these two QTL clusters play major roles. On the other hand, the results of our QTL analysis are not consistent with a previous report that the major QTL for V. anguillarum disease resistance is located on LG11 (LG7 in that study) in Japanese flounder.[41] In the previous study, two V. anguillarum disease resistance-associated polymorphic SSR markers were selected using bulked segregant analysis and then directly anchored onto LG11. The QTLs were finally identified by genotyping of 22 SSRs of LG11 in an F1 population.[41] Because the QTL detection analysis was restricted to a specific LG at the very start, other potential QTLs may have been missed. The other possible explanation for the observed discrepancy is that it was the inevitable result of using a different family under different environmental conditions, just as repeated analysis of the same family treated by two successive challenges would result in lower repeatability for QTL detection for such a complex quantitative trait.[42,43] To enhance the accuracy of QTL detection, results from multiple full-sib families should be integrated in the future.

To identify potential candidate genes, we compared the detected QTLs with the scaffold assembly and annotation of the Japanese flounder reference genome. A total of 12 immune-related genes were discovered, corresponding to the multi-population QTL region. The fact that four of these genes (tap1, satb1, cd40, and cd69) are associated with V. anguillarum disease resistance in Japanese flounder is unsurprising, as they are involved in the molecular function of the MHC class-dependent pathway of antigen presentation. Extensive studies have shown that MHC class I and II molecules play a pivotal role in immune defence systems, because they allow T cells to distinguish self from non-self.[44] In particular, investigations of allelic polymorphism and the pattern of evolution in MHC genes in Japanese flounder have indicated an association between MHC genes and V. anguillarum disease resistance.[28,45-47] These latter results also support the quality of our RAD data, because a key issue for QTL detection and associated gene analysis, given adequate genome coverage, is marker quality.

Comparative genome analysis

Next-generation sequencing technology has enabled generation of draft genomes for individual species as well as rapid development of comparative analyses among multiple genomes based on chromosomal-assembly levels.[48] Nevertheless, the large amount of repetitive DNA sequences and high heterozygosity levels in most species, especially marine species, has hindered assembly accuracy and also prevented direct assembly into chromosomal performance.[49] Consequently, a high-quality genetic map is an appropriate way to anchor scaffolds onto chromosomes. In this study, 24 pseudo-chromosomes of Japanese flounder were constructed by anchoring 246 scaffolds totaling 522.99 Mb, accounting for 95.78% of assembly sequences, using a high-resolution SNP genetic map. Such a high level of chromosomal-level genome map integrity is due to the SNPs, which were abundantly distributed across the genome by using the highest-throughput genotyping technique available.[17] The constructed pseudo-chromosomes of Japanese flounder provide a basic foundation to carry out extensive comparative genomics. In the present study, we performed an in silico analysis of orthologous gene pairs, because genes are evolutionarily conserved relative to intergenic regions.[50] From a general point of view, the orthologues (6,768) of flounder-medaka were more abundant than the orthologues (6,463) of flounder–zebrafish, reflecting the relatively closer phylogenetic relationship of medaka and flounder.[51] Looking more closely at evolutionary relationships based on the chromosomal segments, multi-syntenic correspondences were exhibited between flounder and zebrafish, whereas almost all flounder and medaka chromosomes had a 1:1 correspondence. This difference in synteny directly reflects that the flounder have a phylogenetic closeness to medaka than to zebrafish and also suggests that a much larger number of evolutionarily significant events, such as fusion, breakage, and translocation, have occurred in the zebrafish lineage. The medaka–flounder comparative analysis uncovered little evidence of such inter-chromosomal rearrangements in the medaka lineage despite the difference in chromosome numbers among those fish species (zebrafish has 25 chromosomes, whereas flounder and medaka have 24).[52-54] Furthermore, we inferred the occurrence of many intra-chromosomal events that have interrupted syntenic linkage blocks both between flounder and medaka and between flounder and zebrafish; these interruptions may have been the result of gene loss or divergence after duplication, which can contribute to deviations from synteny.[55] Most teleost species have experienced three rounds of whole-genome duplication. Thousands of duplicated genes remaining after divergence of teleost species have subsequently undergone non-reciprocal loss, pseudogenization, or sub-functionalization, and may be associated with unequal microsynteny in flounder and zebrafish.[56,57] From a different perspective, most chromosomal regions among the three analysed teleost species generally had a conserved syntenic relationship. This observation not only supports the marker ordering on the map but also facilitates the functional inference of genes in flounder and further benefits the analysis of genome evolution and comparative genomics.

In conclusion, we used RAD-seq technology for large-scale identification of SNPs that were then successfully used for high-throughput genotyping and construction of a high-resolution genetic map of Japanese flounder. The developed genetic map is the most comprehensive genetic map to date for this species. Through SNP mapping analysis, we identified 9 positive QTLs for V. anguillarum disease resistance that will be of interest to achieve breeding goals for Japanese flounder. We also anchored the genome scaffolds to pseudo-chromosomes and further identified complex syntenic relationships among flounder, medaka, and zebrafish by comparative genomic analysis. The large numbers of generated SNPs and the dense genetic map, coupled with future re-sequencing of multiple breeding families, should not only lay a foundation for chromosomal-level analysis of the flounder genome but should also provide an excellent resource for future molecular breeding efforts such as genome selection.

Supplementary data

Supplementary data are available at http://www.dnaresearch.oxfordjournals.org.

Funding

The work was supported by grants from State 863 High-Technology R&D Project (2012AA10A408), National Nature Science Foundation of China (31461163005), JSPS RONPAKU (Dissertation PhD) Program for C.S., and Taishan Scholar Project Fund of Shandong of China. Funding to pay the Open Access publication charges for this article was provided by the State 863 HIgh Technology R&D Project (2012AA10A408).