Highly dense linkage maps from 31 full-sibling families of turbot (Scophthalmus maximus) provide insights into recombination patterns and chromosome rearrangements throughout a newly refined genome assembly

Full Paper Highly dense linkage maps from 31 full-sibling families of turbot (Scophthalmus maximus) provide insights into recombination patterns and chromosome rearrangements throughout a newly refined genome assembly F. Maroso 1, *, M. Hermida 2 , A. Milla´n 1 , A. Blanco 2 , M. Saura 3 , A. Ferna´ndez 3 , G. Dalla Rovere 4 , L. Bargelloni 4 , S. Cabaleiro 5 , B. Villanueva 3 , C. Bouza 2 , and P. Martı´nez 2, * 1 GENEAQUA S.L, 2007 Lugo, Spain, 2 Departamento de Zoologı´a, Gene´tica y Antropologı´a Fı´sica, Facultad de Veterinaria, Universidade de Santiago de Compostela, 27002 Lugo, Spain, 3 Departamento de Mejora Gene´tica Animal, Instituto Nacional de Investigacio´n y Tecnologı´a Agraria y Alimentaria (INIA), 28040 Madrid, Spain, 4 Department of Comparative Biomedicine and Food Science, University of Padova, Legnaro, Padova, Italy, and 5 Cluster de Acuicultura de Galicia (Punta do Couso), Agui~ no-Ribeira 15695, Spain *To whom correspondence should be addressed. Tel. þ34 982 822423. Email: [email protected] (F.M.); Tel. þ34982822425. Fax. þ982822001. Email: [email protected] (P.M.) Edited by Dr. Toshihiko Shiroishi Received 15 February 2018; Editorial decision 3 May 2018; Accepted 5 May 2018 Abstract Highly dense linkage maps enable positioning thousands of landmarks useful for anchoring the whole genome and for analysing genome properties. Turbot is the most important cultured flatfish worldwide and breeding programs in the fifth generation of selection are targeted to improve growth rate, obtain disease resistant broodstock and understand sex determination to control sex ratio. Using a Restriction-site Associated DNA approach, we genotyped 18,214 single nucleotide polymorphism in 1,268 turbot individuals from 31 full-sibling families. Individual linkage maps were combined to obtain a male, female and species consensus maps. The turbot consensus map contained 11,845 markers distributed across 22 linkage groups representing a total normalised length of 3,753.9 cM. The turbot genome was anchored to this map, and scaffolds representing 96% of the assembly were ordered and oriented to obtain the expected 22 megascaffolds according to its karyotype. Recombination rate was lower in males, especially around centromeres, and pairwise comparison of 44 individual maps suggested chromosome polymorphism at specific genomic regions. Genome comparison across flatfish provided new evidence on karyotype reorganisations occurring across the evolution of this fish group. Key words: Scophthalmus maximus, 2b-RAD sequencing, linkage mapping, genome assembly, chromosome reorganisations V CThe Author(s) 2018. Published by Oxford University Press on behalf of Kazusa DNA Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permission[email protected] 439 DNA Research, 2018, 25(4), 439–450 doi: 10.1093/dnares/dsy015 Advance Access Publication Date: 11 June 2018 Full Paper Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 1. Introduction The information provided by genomes has become invaluable for research fields that deal with animal production, conservation biology or evolution. 1 High-density (HD) genetic maps provide information on genome organisation by establishing the cartography of thousands of markers, which can aid to disentangle the genetic architecture of productive or evolutionary traits. Moreover, these maps represent an essential tool for comparative genomics and can facilitate whole genome assembly. 2,3 The advent of next generation sequencing (NGS) has dramatically reduced the cost and time required for genomic analysis, 4,5 thus contributing to a fast increase of genomic resources. Nevertheless, while DNA sequencing has become an automated task at platforms, genome assembling still relies on the hierarchical ordering of short DNA fragments (contigs, scaffolds), currently being facilitated by the novel long-read assembling approaches (Nanopore, PacBio, Dovetail, Bionano). 6 HD genetic mapping can assist the process of genome assembly through positioning and ordering scaffolds or validating bioinformatic assembly. 7,8,9 This is a particularly challenging task in teleosts, considering the ancestral teleost-specific whole genome duplication event in this group of fish. 10 A huge number of markers and samples can now be analysed at very low cost using genotyping by sequencing (GBS) methods, 4,5,9 which enables to simultaneously identify and genotype thousands of single nucleotide polymorphisms (SNPs) across samples. A popular GBS method is Restriction Associated DNA sequencing (RADseq), a technique that combines the power of NGS with the simplification of genomes through restriction enzyme digestion. 11–13 RADseq is being increasingly used for the construction of HD genetic maps in fish comprising thousands of markers (e.g. Refs 14–19 ). These maps have allowed the identification of genomic regions associated with productive traits 20,21 or with adaptation to specific environments. 14,16 The pattern of samples and families used in fish for map construction has been quite similar, and only very recently, the number of families has increased. 19,22 However, when many families have been analysed, the number of offspring per family has been much smaller (e.g. Tsai et al. 2016 23 : 60 families with less than 10 offspring each; Palaiokostas et al. 2016 24 : 75 families with 10 offspring on average). This is a way to achieve a consensus reference map at species level by averaging inter-individual variation, but the low family sampling size can result in suboptimal estimation of inter-marker distances and mapping order. 25,26 The presence of genomic reorganisations and variation of mapping parameters across individual or sex-specific maps has been hardly addressed to date in fish despite its importance from both evolutionary and practical perspectives. For example, inversions might aid to maintain blocks of coadapted genetic variants at specific environments, such as reported in fish like cod 27 or stickleback. 7,28 Moreover, sex-specific mapping features are important when tracking markers associated to economically relevant traits since recombination patterns in males and females may differ and thus, they should be considered for accurate evaluations. NGS technologies now enable studying intra-specific variation on genome organisation and recombination with enough depth to improve the analysis of associations with traits of interest. This would aid for a more accurate mining to identify responsible genes and mutations and choose the best strategy for marker assisted selection. The turbot (Scophthalmus maximus) is the main worldwide aquaculture flatfish, 9 People’s Republic of China being currently the most important producer (>60,000 tons/year 29 ). Genetic breeding programs in this species began in the 90s. Nowadays, a broad battery of genetic markers and genomic tools are available in turbot, which can aid to understand the genetic basis of the main economically important traits (growth rate, disease resistance, sex control) and for applying a more efficient selection. 30–35 The first genetic maps in turbot, based on anonymous or gene-associated microsatellites and SNPs, 36–38 were later integrated in a consensus map. 38 More recently, a high-density map was reported by Wang et al., 21 but the scarce information provided on genetic marker-associated sequences made it not feasible to integrate with previously reported maps. 9 Centromeres have been located in the turbot genetic map using half-tetrad analysis 40 , and recently, a high coverage draft genome has been assembled and anchored to the reported consensus map. 36 Moreover, the development of genomic resources for other important flatfish species 41–43 and their integration in a common framework has allowed comparative analyses to transfer genomic information on production-related traits between species and to look for insights into flatfish evolution. 9,38 In this study, we genotyped 18,214 SNPs using the 2b-RAD technology in 1,268 individuals from 31 turbot full-sibling (sib) families with the aim of: (i) obtaining a robust consensus genetic map for the species; (ii) integrating the genetic and physical maps in a common framework to facilitate transference of segregation data into the physical genome; (iii) refining the turbot genome assembly and comparative mapping with other flatfish; and (iv) assessing the variation in genome organisation and recombination between sexes and individuals. Our results rendered much deeper information on intra-spe- cific mapping variation than previously reported in fish and provided useful integrated genomic resources to understand the architecture of economically important traits for further applications in industrial breeding programs. 2. Materials and methods 2.1. Families for mapping The broodstock of CETGA (Aquaculture Cluster of Galicia, Ribeira, Spain), a population representative of turbot from the Atlantic area, was used to found a set of 44 families within the goals of the FISHBOOST project (EU 613611). This project is focused on improving the efficiency and profitability of finfish European aquaculture by fostering selective breeding to a next level, including the application of genomic selection. This goal also involved the refinement and improvement of existing mapping resources and their integration in a common framework within the reported turbot genome assembly. 34 The experiment was performed under approval of the ethics committee of CETGA to adjust to animal welfare legislation. A total of 31 full-sib families were selected among the 44 families considering both successful genotyping and informativeness (Fig. 1). The partly factorial mating of 22 males 22 females performed made available 27, either maternal of paternal, half-sibs families (eight males and nine females contributed to more than one family). A total of 1,224 offspring were used for linkage analysis representing on average 39.5 individuals per family (range: 36–45). 2.2. SNP calling and genotyping 2.2.1. DNA extraction and library preparation Genomic DNA was extracted from fin clips using SSTNE buffer (a TNE buffer modified by adding spermidine and spermine) and a standard NaCl isopropanol precipitation. 44 Library preparation followed the 2b-RAD protocol 45 slightly modified to include replicated 440 Highly dense turbot population mapping Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 polymerase chain reactions (PCRs) per sample (to reduce PCR duplicate bias) and to quantify DNA sampling before pooling (to increase sample coverage homogeneity). Detailed information about the protocol can be found in Supplementary Material text. Briefly, DNA of each sample was separately digested using AlfI restriction enzyme, adaptors including specific sample barcodes ligated, and the resulting fragments amplified. After PCR purification, samples were merged in pools of 128 samples for sequencing and genotyping the offspring, while the 44 parents were genotyped in a single pool, to increase coverage and genotyping accuracy. A total of 11 pools (1 for the parents and 10 for the offspring) were sequenced on an Illumina NextSeq sequencer using a 50 bp single end protocol at BMR Genomics (Padova, Italy). 2.2.2. Raw data filtering and bioinformatics Raw data were demultiplexed at the sequencing facility. All reads were trimmed at 36 bp and centred on the enzyme recognition site. Stacks’ process_radtag pipeline 46,47 was used to filter raw reads according to sequencing quality. All reads with an uncalled base were removed, and quality was checked using a sliding window of 25% of the sequence length. If the mean quality score of this window dropped below 20, the sequence was removed. Filtered reads were mapped against the turbot genome 34 using Bowtie 1.2.0 48 and allowing a maximum of three mismatches in the first 33 bases. Only a unique match on turbot genome was permitted, otherwise reads were discarded. Mapped reads from parent samples were fed into Stacks’ module ref_map.pl using the bounded model with a significance threshold of 0.05 and an upper bound of 0.05. As a further filter, only SNPs genotyped in 80% of the parents with a minimum read depth of 10 reads were retained using the populations module. When a RAD-tag contained more than one SNP, only the first SNP from the 5’ end was retained. The set of SNPs obtained from parents was used as a reference for offspring genotyping using the Stacks module sstacks. The obtained dataset (parents and offspring) was further filtered according to population parameters and Mendelian segregation using the 26 most informative families. Those markers showing consistent deviations from Mendelian segregation (P<0.01) across families or extreme departures from Hardy– Weinberg equilibrium (HWE) in the parental population (P<0.001) due to heterozygote excess (F IS <0.5, P<0.0001; paralogous genes) or deficit (F IS >0.5, P<0.0001; null alleles) were filtered out. Finally, RAD-tags containing SNPs were mapped against the turbot gene catalog 38 using Blastn (E-value <e20) to obtain sequence annotations. 2.3. Construction of genetic maps Genotypes of filtered SNP markers were properly coded as JOINMAP Cross Pollinator (CP) cross type with unknown linkage phase. Markers genotyped in less than 90% of the offspring were excluded for linkage analysis in each family. All genetic maps (individual, family, male consensus, female consensus and species consensus) were constructed using JoinMap 4.1. 49 First, markers were associated to their linkage groups (LGs) using the Grouping function of JoinMap based on a series of LOD scores (logarithm of the likelihood ratio that two genes are linked regarding that they segregate independently) increasing by one, from 3 to 10. The LOD score was then selected in each family based on the number of turbot chromosomes (n¼22); in all cases, a LOD >6.0 was used. Unlinked markers and markers assigned to very small LGs exceeding the haploid number of 22 chromosomes were excluded from further analysis. Second, marker ordering was performed using the Maximum Likelihood algorithm with default parameters, and the recombination frequencies converted into map distances in centi-Morgans (cM) using the Kosambi mapping function. For each family, individual Figure 1. Families selected for linkage analysis in S. maximus. Male (M) female (F) crosses and the number of offspring per family are indicated. Note the presence of several half-sib families sharing the father or the mother. 441F. Maroso et al. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 female and male maps, as well as a consensus family map, were generated. Due to computational limitations, only the three most informative families were selected to construct the species consensus, the female consensus and the male consensus maps using MergeMap. 50 All maps constructed were drawn using MapChart ver. 2.3. 51 In order to investigate the heterogeneity of recombination frequency (RF) between sexes, RFs were extracted from JoinMap analysis. The pattern of RF in the father and the mother could not be straightforwardly compared within each family due to the different panel of informative SNPs; consequently, the set of markers shared between both parents was selected in each family for a detailed evaluation of RF variation between sexes across LGs within each family (heterozygous in both parents with 1:2:1 segregation). 2.4. Genetic versus physical map: anchoring scaffolds to linkage groups The consensus genetic map was used to refine the turbot genome assembly 34 and to obtain a physical map consistent on 22 megascaffolds, representing the complete sequence of each of the 22 LGs of turbot map. 32,39 This task implied (i) to assign previously unlinked scaffolds to specific LGs in the new high-dense (HD) genetic map; (ii) to correct scaffold misassemblies by comparing the genetic and physical maps; and (iii) to establish the correct orientation of scaffolds within each megascaffold. All contigs including at least one marker significantly assigned to a LG were anchored in the new map. When markers of a specific scaffold were assigned to more than one LG, the scaffold was split into fragments pertaining to different LGs or to the same LG but in different positions (see Results). To establish the most confident splitting of these scaffolds, we compared the sequence of the gap between the flanking markers assigned to different LGs with the two model fish genomes, Gasterosteus aculeatus and Oryzias latipes, which had demonstrated large synteny with turbot. 34,38 Then, we looked for the turbot genes in the gap assigned to different LGs in the model species in order to narrow down as much as possible the gap (Supplementary Fig. S1). Following this approach, the narrowest gap was finally cut in the middle and halves assigned to each LG. Reorganisations were only validated when information was consistent with both model species. To find out the correct order and orientation of scaffolds within LGs, we obtained the average map position of the markers belonging to each scaffold using the maps of the 10 most informative families. Then, we checked for the correlation between the physical position of markers and their genetic position. When the correlation was positive, the original orientation was maintained, and when negative, the reverse complementary sequence of the scaffold was retrieved for assembling. 2.5. Integrating genetic and physical maps with cytogenetic data Genetic and physical maps were then integrated with the turbot cytogenetic map using previous information. 32,40,52 Starting from centromere positions 40 and considering the relationship between genetic and cytogenetic maps, 32 we established the correspondence between LGs, megascaffolds and chromosomes as far as possible. Then, we compared the consensus genetic map of each chromosome with its physical assembly (megascaffold) to assess their correspondence with centromere, telomere and interstitial heterochromatin locations. When no marker was available close to the centromeric region, the closest markers around the centromere were used as reference or, in a few cases, the farthest marker from the centromere was used for its rough location at the other chromosome end, considering the acrocentric morphology of most turbot chromosomes. 52 Finally, this comparison was also carried out between sexes using the sex-specific consensus maps, in order to assess RF differences across chromosomes between males and females. 2.6. Inter-individual genomic reorganisations in turbot The availability of 31 family maps and 44 parental maps (Fig. 1), enabled an exhaustive comparison of genetic maps between individuals and sexes to assess inter-individual variation on RF across the genome. This variation, which might be related to chromosome polymorphisms existing in turbot populations, had not been previously addressed in fish with the number of markers and families here evaluated. In order to identify putative genomic reorganisation between individuals, we performed a pair-wise comparison between all parental maps by plotting the positions of each marker. Plots were created with R 53 and the different LGs were labelled with different colours for a better visualisation. In order to handle the amount of data obtained, 62 individual maps at 22 LGs (due to half-sib families; Fig. 1), each LG was divided in stretches of 50 SNPs representing 1.6 Mb on average (e.g. 23 blocks for LG02, the longest LG with 1,132 SNPs). Then, the correlation between the physical and genetic maps was analysed for each block at each individual map in order to identify potential inversions. A heat map was generated to provide a more straightforward view of the results. 2.7. Interspecific genomic reorganisations between turbot and other flatfish species The improved turbot genome was further exploited to complete and refine previous information about genomic reorganisations and chromosome evolution of flatfish. 9,34 For this, we compared the improved turbot genome with those of Japanese flounder (Paralichthys olivaceus 43 ) and half-smooth tongue sole (Cynoglossus semilaevis 42 ) at chromosome level using LASTZ 54 (using options ‘–notransition – step ¼400 – nogapped – format ¼rdotplot’) and compared each turbot chromosome with the orthologues of the other two species. Further, in order to assess the relation between flatfish megascaffolds to a model fish genome outside Pleuronectiformes, an additional LASTZ analysis was carried out with stickleback (Gasterosteus aculeatus) using the same parameters. 3. Results 3.1. SNP calling, filtering and genotyping A total of 3,946,354,960 reads (3.9 Gb) were produced at the BMR sequencing platform, corresponding to 1,282 individuals represented by 3,078,280 reads on average (range: 312 – 9,610,380). Fourteen offspring were removed due to the low number of usable reads (<75,000), leaving a total of 1,268 individuals (44 parents and 1,224 offspring) pertaining to 31 full-sib families for subsequent analyses. After quality filtering and alignment to the turbot genome, the average number of reads per sample was 6,663,809 (range: 2,601,955– 8,946,621) for parents and 2,774,579 (range: 283,203–5,596,610) for offspring. After Stacks’ populations filter, 139,293 RAD-tags were retained, of which 25,511 (18.3%) contained at least one SNP. Since only one SNP was retained per RAD-tag, a total of 25,511 SNPs were considered after eliminating the remaining SNPs in the same RAD-tag (Supplementary Table S1). Among these 25,511 442 Highly dense turbot population mapping Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 SNPs, 18,214 were finally retained after filtering loci by: (i) consistent deviations from Mendelian segregation (5,866 loci: P<0.01 in at least two families and representing more than 20% of the informative families); (ii) very low minimum allele frequency in the parental population (<0.015) suggestive of genotyping errors (991 SNPs); and (iii) extreme deviations from HWE suggesting either paralogous genes (364 loci) or high null allele frequency (91 loci) (Supplementary Table S1). This set of 18,214 SNPs represented a frequency of one SNP every 163 bp in the reduced portion of the genome sequenced with 2b-RAD, a frequency slightly lower than that reported in other studies in turbot, 55,56 and a density of one SNP every 31,130 bp in the whole turbot genome (567 Mb 34 ). The average proportion of genotyping errors (Mendelian inconsistencies in families) of our 2b-RAD genotyping approach was 0.68% throughout all families in the filtered SNP dataset. Out of the 18,214 SNPs finally retained, nearly half matched to turbot genes (9,028). Among those lying within exons (834), 400 corresponded to non-synonymous and 430 to synonymous allelic variants (Supplementary Table S1). 3.2. Construction of genetic maps On average, 5,385 SNPs were used for map construction in each of the 31 families analysed (range: 2,079 in Fam55 to 6,757 in Fam30; Supplementary Table S2). Maps consisted of 22 LGs corresponding to the 22 chromosomes of the turbot haploid karyotype, averaging 140.5 cM and 245 markers per LG and family (Supplementary Tables S2 and S3). Note that LGs nomenclature follow previous publications 32 and accordingly, LG18 is missing due to its merging with LG08 (LG08 þ18, here LG08). LG22 was on average the shortest LG (94.9 cM; range: 28.1 cM in Fam55 to 223.5 in Fam19) and included the lowest number of markers (189.9; range: 56 in Fam16 to 257 in Fam39). LG02 was on average the longest LG (207.7 cM; range: 49.9 cM in Fam59 to 414.0 cM in Fam19) with the highest number of markers (341.7; range: 97 in Fam55 to 457 in Fam39). A total of 17,169 SNPs were successfully mapped at least in one family, 1,045 remaining unmapped. Of the last, 1,029 were located in scaffolds where other markers were anchored to the linkage map and could therefore be assigned to specific LGs. Only 16 SNPs could not be assigned to any LG. Marker interval across all families averaged 0.6 cM and ranged from 0.43 cM in Fam11 to 1.08 cM in Fam19. If unique positions instead of markers were accounted for, marker interval averaged 2.0 cM and ranged from 1.81 in Fam13 cM to 2.74 cM in Fam55. The turbot consensus map constructed with the three most informative families using MergeMap included 11,845 markers distributed over 22 LGs representing a total length of 8,532.6 cM (Supplementary Table S3A;Fig. 2). LG length varied from 282.9 cM (LG22) to 588.7 cM (LG02), with an average length of 387.8 cM. The average marker interval was 0.72 cM and ranged from 0.57 cM in LG01 to 0.89 cM in LG4. Because of the inflated length observed for this consensus map, we also estimated the average LG lengths from the three families used to construct the consensus, as these values represent more realistic lengths than those rendered by the consensus MergeMap algorithm (see Discussion) (Table 1). These LG lengths sum up to 3,753.9 (i.e. an average marker interval of 0.32 cM). Female and male consensus genetic maps differed significantly totalling 8,274.8 cM and 6,791.6 cM (F:M ratio ¼1.22), with average LG length of 376.1 and 308.7 cM, respectively. In both females and males LG02 was the largest LG (545.8 vs 474.9 cM) and LG22 the shortest one (292.1 vs 198.7 cM) (Supplementary Fig. S2). As for the species consensus map, Table 1 reports average LG lengths also for female and male consensus maps. Individual female and male genetic maps (62 in total; Supplementary Tables S2 and S3) were also constructed for a more detailed comparison of RF between individuals and sexes, and to look for putative chromosome polymorphisms in the turbot population (44 parents). The length of female genetic maps was higher than that of male maps in all families and for all LGs (average F:M ratio ¼1.36). For a more detailed comparison, 10,168 markers shared between sexes and mapped in at least one family were selected (Supplementary Table S4). As expected, the average inter-marker distance for the same pair of markers was significantly higher in females than in males (2.2 and 1.7 cM, respectively; P<0.001, t-test for related samples). When F:M intervals were compared at each LG, the difference was significant for all groups excluding LG02, LG17 and LG22 (P>0.05) (Supplementary Table S4). When F:M intervals were compared across all LGs within each family, again inter-marker distances were higher and significant in females in all but four families (families 16, 22, 47 and 55, P>0.05; Supplementary Table S4). 3.3. Anchoring turbot genetic and physical maps: megascaffolds Scaffolds were assigned to the 22 LGs of the turbot consensus map, then ordered and oriented, thus constituting 22 megascaffolds. Series of 100 Ns were used in the megascaffold sequence to show the gap between adjacent scaffolds denoting missing information. Additionally, the short uncertain regions around the boundaries of split scaffolds were written in lowercase to highlight their lower confidence. The complete nucleotide sequence of the 22 megascaffolds with their correspondent genes located and annotated has been uploaded to NCBI (SUB3239805). A total of 28 scaffolds were split and reallocated according to the very consistent mapping information in all the 31 families analysed (Supplementary Table S5). Most scaffolds (22) were split into two fragments and reallocated to different LGs. In three cases, scaffolds were split in three different LGs (sm5_s00043, sm5_s00056 and sm5_s00084), and in other three cases, the two fragments of the scaffold were reallocated within the same LG (sm5_s00006, sm5_s00007 and sm5_s00049). In all cases, the boundaries of the breakage point were assigned with notable accuracy (Supplementary Fig. S1), and uncertain gap regions averaged 38,010 bp, ranging from 549 bp (sm5_s00086) to 160,155 bp (sm5_s00084) (Supplementary Table S5). Rearranged scaffolds represented a total of 108.5 Mb, to say 20% of the entire turbot genome. A total of 17,169 SNPs corresponding to 558 scaffolds (after splitting 28) rendered a significant match to the turbot genetic map (Supplementary Table S1). On average each scaffold contained 32.7 SNPs ranging from 610 SNPs in sm5_s00003_17 (LG17) to one SNP in 204 scaffolds. The total length of the anchored sequences was 524,925,810 bp, representing 96.5% of the turbot genome assembly (544 Mb). The average megascaffold length was 23,860,264bp, ranging from 14,917,408 bp (LG22) to 31,897,078bp (LG02) (Table 1). Megascaffolds contained on average 827 SNPs (range 578–1144) and 25.4 scaffolds (range 14–46) (Table 1). All scaffolds assigned to chromosomes were consistently ordered within the megascaffold according to the mapping information across families; the orientation within each chromosome could be verified for those 336 that contained at least two informative markers. The remaining 204 scaffolds with only one marker could not be oriented, but they 443F. Maroso et al. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 represented a very minor proportion of the total assembly (2.0%) because of their small size. Interestingly, they were mostly located at the end of chromosomes (telomeres), suggesting assembling problems caused by repetitive motifs. 3.4. Integrating cytogenetic information: chromosomes vs LGs vs megascaffolds Turbot genetic, physical and cytogenetic information obtained in this work and from previous studies was integrated as far as possible, and accordingly, a unified nomenclature is proposed in this study (Table 1). Four chromosomes (C1, C2, C3 and C22) 34 were unequivocally assigned to specific LGs: the two metacentric C1 (LG02) and C2 (LG16), the submetacentric nucleolus organiser region (NOR)-bearing chromosome C3 (LG15) and the smallest acrocentric chromosome C22 (LG22), constituting the four marker chromosomes identifiable in the turbot karyotype. The remaining megascaffolds, which correspond to acrocentric chromosomes of decreasing size, were termed according to their physical length from C4 to C21 (Table 1). The consensus genetic and physical maps showed a linear correspondence along most chromosome length in all LGs, excluding those regions around centromeres and telomeres (Fig. 3A and Supplementary Fig. S3). A cloud of dots indicative of poor mapping accuracy were observed around most mapped centromeres 40 (i.e. LG01, LG10, LG11, LG19 and LG20). Dot clouds were located in intermediate regions of the metacentric chromosomes (LG02 and LG16) or at the terminal regions in the acrocentrics in accordance with their cytogenetic structure. LG15, the NOR-bearing chromosome, showed a moderate dispersion of dots throughout the proximal centromere chromosome half likely related to the tandem rDNA clusters and their associated heterochromatin in that region. 52 Also, close to telomeres, although not so markedly, dots were either Figure 2. Consensus genetic map of S. maximus using the three most informative full-sib families; genetic distance in cM in the left bar; the position of centromeres (C) or genetic markers far from centromere (T) in acrocentric chromosomes are indicated. 444 Highly dense turbot population mapping Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 dispersed (i.e. LG03, LG08, LG11 and LG14) or showing particular drawings (LG02, LG06 and LG15) suggestive of wrong physical assembly or genetic mapping. Finally, the pattern of the acrocentric chromosome C5 (LG04) was remarkable because of the presence of an internal cloud of dots in addition to that associated around the centromere at the chromosome end. This LG very likely correspond to the chromosome 5 of turbot karyotype, which showed an interstitial heterochromatic block with C-banding. 52 When the same comparison between genetic and physical maps was performed for the female and male consensus maps, dot clouds around centromeres were very prominent in the male map and subtle, and sometimes inexistent, in the female maps (Fig. 3B and Supplementary Fig. S4; i.e. LG01, LG05, LG09, LG15, LG19 and LG22). So, it seems that the length difference between male and female maps is mainly related to the lower recombination rates around centromeres characterising male maps. 3.5. Intraspecific chromosome reorganisations Marker positions showed a good correspondence between the 62 parental genetic maps (including ‘replicates’ from half-sib families) although particular drawings were suggestive of different marker Table 1. Average length and constitution (SNP and scaffold number) of LGs in the species, female and male S. maximus consensus maps Karyotype Physical map Species consensus map Female consensus map Male consensus map Chromosome ID LG Megascaffold length (bp) No SNPs No scaffolds Average contig length (bp) Mapped markers Average length Mb/cM Mapped markers Average length Mb/cM Mapped markers Average length Mb/cM C7 01 26,874,605 1,031 32 839,831.4 671 178.3 0.151 487 176.6 0.152 409 162.9 0.165 C1 02 31,897,078 1,144 46 708,824.0 761 255.8 0.125 513 261.6 0.122 546 215.0 0.148 C17 03 21,400,935 785 23 930,475.4 506 180.7 0.118 317 179.1 0.119 395 156.3 0.137 C5 04 29,497,696 907 31 951,538.6 610 220.8 0.134 412 251.0 0.118 439 171.7 0.172 C12 05 24,809,483 707 20 1,240,474.2 443 125.0 0.198 300 129.1 0.192 275 107.2 0.231 C10 06 25,205,370 851 26 969,437.3 587 186.1 0.135 395 205.3 0.123 455 153.3 0.164 C13 07 24,510,356 699 20 1,225,517.7 412 168.3 0.146 317 179.4 0.137 301 143.3 0.171 C4 08 30,930,439 1,066 38 813,958.9 640 266.3 0.116 434 303.3 0.102 427 196.0 0.158 C8 09 25,788,341 929 20 1,289,417.1 593 200.0 0.129 415 209.6 0.123 406 165.1 0.156 C11 10 25,098,921 922 25 1,003,956.8 650 169.1 0.148 419 166.8 0.150 497 156.5 0.160 C6 11 27,408,095 808 33 830,548.3 494 150.9 0.182 406 174.1 0.157 323 108.9 0.252 C9 12 25,240,469 796 21 1,201,927.1 536 194.9 0.130 395 204.2 0.124 376 162.2 0.156 C19 13 20,256,549 726 20 1,012,827.5 472 154.0 0.132 349 176.3 0.115 346 119.2 0.170 C16 14 21,450,603 778 20 1,072,530.2 514 145.1 0.148 374 158.2 0.136 402 114.3 0.188 C3 15 24,127,324 887 33 731,131.0 593 146.4 0.165 413 156.7 0.154 383 113.0 0.214 C2 16 23,833,596 813 32 768,825.7 531 197.8 0.120 375 201.1 0.119 371 171.1 0.139 C21 17 15,881,475 680 16 992,591.8 467 127.6 0.124 335 141.3 0.112 332 103.9 0.153 C15 19 21,781,057 755 22 990,048.0 492 128.7 0.169 372 123.0 0.177 345 118.0 0.185 C14 20 22,751,219 776 23 989,183.4 456 145.5 0.156 304 163.1 0.139 322 116.4 0.195 C18 21 21,354,662 760 21 1,016,888.7 501 136.2 0.157 356 138.6 0.154 386 123.8 0.172 C22 22 14,917,408 578 14 1,065,529.1 408 121.0 0.123 291 128.6 0.116 259 102.4 0.146 C20 23 19,910,129 796 22 905,005.9 508 155.4 0.128 378 151.4 0.132 378 129.4 0.154 Sub-total 524,925,810 18,194 558 11,845 3,753.9 0.140 8,357 3,978.4 0.132 8,373 3,109.9 0.169 Not anchored 19,312,992 16 15,936 1,211.9 6,369 9,857 9,841 Total 544,238,802 18,214 16,494 18,214 18,214 18,214 Figure 3. Correspondence between physical (abscissae) and genetic (ordinates) maps showing the dispersion of dots around the centromere at LG11 (A) and the difference between male and female maps at LG11 (B). The vertical bar in both figures show the position of the centromere according to Mart ınez et al. (2008) 42 . Genetic markers corresponding to different scaffolds are shown in different colors. 445F. Maroso et al. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 order between individual maps (Fig. 4 and Supplementary Fig. S5). The negative correlations between physical and genetic maps at these regions supported this observation, and the presence of putative chromosome rearrangements that involved either single blocks (50 SNPs, 1.6 Mb) or longer chromosome tracts was identified (Supplementary Fig. S6). Inverted single blocks were very frequent at LG02 (blocks 7 and 23), LG05 (block 12), LG08 (block 2), LG13 (blocks 14 and 15) and LG22 (block 12), while larger inverted chromosome tracts were observed in some individuals at LG06 (blocks 5–9) and LG09 (blocks 13–16) (Supplementary Figs S5 and S6). However, these putative inversions, particularly those related to single blocks, were found close to telomeres or centromeres, where mapping accuracy is poorer both due to assembling problems and to the lower ratio of genetic/physical distance. This suggests that they are more likely related to technical devices than to true reorganisations. Otherwise, the longer tracts observed at LG06 and LG09 were located at interstitial regions, and thus, signatures of inversions should be more consistent. 3.6. Interspecific chromosome reorganisations The whole turbot megascaffold sequences were aligned against the genomes of Japanese flounder and tongue sole using the LASTZ program and a 1:1 macrosyntenic pattern was observed (Supplementary Fig. S7A and B), as reported previously. 9 However, some Robertsonian translocations (fusions or fissions) and other minor inter-chromosomal rearrangements were observed. The most notable cases when comparing turbot (n¼22) and Japanese flounder (n¼24) genomes were the correspondence of turbot LG02 with Japanese flounder Chr06 and Chr14 (Fig. 5A) and that of LG16 with Chr09 and Chr16 (Supplementary Fig. S7A). In the same way, turbot LG02 was related to tongue sole Chr13 and Chr14 (Fig. 5B), which suggests that LG02 fusion is a derived condition (apomorphy) in the turbot lineage. The position of turbot centromeres at LG02 was the expected one if a centric fusion between two chromosomes had occurred in the turbot lineage (Fig. 5), as also occurred for LG16 (Supplementary Fig. S7A). These chromosome fusions, starting from the ancestral teleost karyotype (2n¼48), 57 would explain the karyotype differences between turbot and Japanese flounder, as suggested by Robledo et al. 9 Despite the general macrosyntenic pattern observed, pairwise alignment of chromosomal sequences indicated that extensive intrachromosomal reorganisations (mainly inversions) would have taken place throughout the evolution of these species (Supplementary Fig. S7A and B). Reorganisations were much more extensive when comparing turbot to half-smooth tongue sole than to Japanese flounder (Fig. 5C and D). The results of the comparison with stickleback are presented in Supplementary Figure S7C. Similar fusion patterns as those identified in the comparison with the other flatfish species are visible, as for example at LG16, that showed macrosyntenic patterns with stickleback’s ChrVII and ChrIV. 4. Discussion The reduced cost of the new genotyping by sequencing (GBS) techniques has facilitated the development of highly dense genetic maps in aquaculture fish species, 3,17–19 including turbot. 21 Here, we combined for the first time a high number of markers and families for a deep mapping analysis aimed at integrating genetic and physical map resources in turbot to understand its genome organisation and variation, a study not performed to date in fish. 23,58 The average number of SNPs in the 31 full-sib families analysed in our work was in the upper range reported in the most recent linkage mapping studies (5–6k SNPs) and 40 individuals were genotyped at each family. The turbot consensus map (8,532.6 cM) was much longer than any previous linkage map reported for the species (6,647 SNPs, 2,622.1 cM 21 ), but this elongation is the consequence of combining a high number of markers with several mapping families and the MergeMaps software, an artefact highlighted previously by different authors. 24,59–61 A direct consequence of this ‘inflated’ map length is that the average ratio between physical and genetic distance is much lower than usual (0.066 Mb/cM vs 0.33 Mb/cM in teleost). A way to compensate this bias is to normalise for each LG consensus map coordinates to the average genetic distance obtained from the individual maps used to construct the consensus 62 (see Table 1). Accordingly, the average length and average inter-marker distance of the three single family maps used to construct the consensus was 3753.9 and 0.32, respectively, figures much closer to those found by Wang et al. 21 using RADseq. Interestingly, the sex-determining region-bearing chromosome (LG05 30 ) showed one of the lowest RF ratio (cM/Mb), which agrees with theoretical expectations on the evolution of the sex-determining chromosome. 32 The anchorage of the assembled turbot genome to a robust and dense genetic map is essential to transfer marker associations detected by segregation analysis in families to physical positions in the genome for gene mining. 63 In our study, 20% of turbot genome assembly 34 was refined, thanks to the mapping information obtained across the 10 most informative families. Furthermore, scaffolds representing 96.5% of turbot genome assembly were anchored to chromosomes, ordered and oriented, and thus, a megascaffold was set up for each LG. This represents a significant improvement from the 80% anchored with the previous medium density turbot map. 34 Also, the available cytogenetic information 32,40,52 was integrated as much as possible in a common framework with genetic and physical maps, and four marker chromosomes of the turbot karyotype were associated with their respective megascaffolds, and the remaining 18 roughly assigned according to their physical length. Using this information, a new nomenclature is here proposed for the 22 chromosomes/LGs of the species, which will facilitate the integration of the different levels of genome organisation for comparative studies within flatfish and teleosts. 9 The high correlation observed between genetic and physical maps supports the confidence of marker ordering and physical assembling of turbot genome. However, this general trend was lost at specific chromosome structures such as centromeres and telomeres. Around Figure 4. Particular drawings denoting the lack of collinearity between two individual genetic maps of S. maximus: Female Fam09 vs Male Fam07 maps; consecutive LGs are represented in different colours from LG01 (left bottom) to LG23 (right top). 446 Highly dense turbot population mapping Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020 centromeres extensive dot clouds were visualised, suggestive of poor mapping accuracy very likely related to the decay of RF in the vicinity of heterochromatic blocks. 64,65 This was observed not only for pericentromeric heterochromatin but also for interstitial heterochromatic blocks such as the most conspicuous one detected at LG04 (C5). A remarkable RF difference was observed between males and females around heterochromatic structures, especially the centromeres, where females hardly showed RF decay. This fact is the main responsible for the length difference between male and female maps of turbot, as reported previously in other fish species. 58,66–68 It is tempting to ascribe the different RF between males and females observed in fish 19,66 to different crossing-over patterns around heterochromatic blocks. However, other explanations have been proposed for different species, such as the GC content in Asian sea bass 19 ; stronger selection pressure in male gametes during the haploid life stage in zebrafish 66 or in some plants 69 ; and differences in chromatin distribution between sexes as for example in humans 70 . More information is needed for a comprehensive explanation to these observations, and indeed, NGS technologies can provide suitable data for this task. The lack of correspondence between physical and genetic maps observed at telomeres in several families and LGs is likely a consequence of a wrong assembly and/or mapping at chromosome ends. These regions were characterised by the presence of a number of short contigs, suggestive of assembling problems likely due to the presence of repetitive DNA. Similarly, when many markers and a Figure 5. LASTZ plots between S. maximus linkage groups (LG) and chromosomes (Chr) of P. olivaceus and C. semilaevis: (A) LG02 vs Chr06 and 14 of P. olivaceus; (B) LG02 vs Chr13 and 14 of C. semilaevis; (C) LG04 vs Chr07 of P. olivaceus; (D) LG04 vs Chr05 C. semilaevis. The gross vertical bar points the position of centromeres according to Mart ınez et al. (2008) 42 and the light ones the limits between consecutive scaffolds. 447F. Maroso et al. Downloaded from https://academic.oup.com/dnaresearch/article-abstract/25/4/439/5035912 by guest on 29 May 2020