Draft genomes of two Atlantic bay scallop subspecies Argopecten irradians irradians and A. i. concentricus

The two subspecies of Atlantic bay scallop (Argopecten irradians), A. i. irradians and A. i. concentricus, are economically important aquacultural species in northern and southern China. Here, we performed the whole-genome sequencing, assembly, and gene annotation and produced draft genomes for both subspecies. In total, 253.17 and 272.97 gigabases (Gb) of raw reads were generated from Illumina Hiseq and PacBio platforms for A. i. irradians and A. i. concentricus, respectively. Draft genomes of 835.7 Mb and 874.82 Mb were assembled for the two subspecies, accounting for 83.9% and 89.79% of the estimated sizes of their corresponding genomes, respectively. The contig N50 and scaffold N50 were 78.54 kb and 1.53 Mb for the A. i. irradians genome, and those for the A. i. concentricus genome were 63.73 kb and 1.25 Mb. Moreover, 26,777 and 25,979 protein-coding genes were predicted for A. i. irradians and A. i. concentricus, respectively. These valuable genome assemblies lay a solid foundation for future theoretical studies and provide guidance for practical scallop breeding.


Background & Summary
Two subspecies of the Atlantic bay scallop, the northern subspecies Argopecten irradians irradians (Lamarck, 1819) and the southern subspecies A. i. concentricus (Say, 1822), are widely cultured in China. The northern subspecies A. i. irradians is mainly cultured in northern waters, while the southern subspecies A. i. concentricus is generally cultured in southern waters 1 . Both subspecies were introduced from the USA between the 1980s and 1990s. In general, these bay scallops grow fast but have short life spans (i.e. <24 months) 2 . These two subspecies are morphologically similar, although the ratio of shell width (W) to shell height (H) or shell length (L) of A. i. concentricus is remarkably higher than that of A. i. irradians-the average W/L ratio of adult A. i. concentricus and A. i. irradians are 0.59 and 0.45, respectively 3 . Although both subspecies are adapted to their natural habitats with wide temperature ranges, A. i. irradians is more tolerant to the cold northern waters but cannot survive in the southern warm waters, whereas A. i. concentricus is better adapted to the warm southern waters but stops growing at a temperature of 12 °C or lower.
Successful diallel crossbreeding has been performed between the two subspecies, as well as between the Peruvian scallop (Argopecten purpuratus) and both of the two bay scallop subspecies 4,5 . In addition to high fertilization and hatching rates, the resulting F 1 hybrids exhibited excellent performance in production traits such as growth and survival, indicating a great potential in stock improvement via inter-or intra-specific hybridization between different subspecies or populations. To date, three high-performance strains, 'Bohai Red' , 'QN-2' and 'QN Orange' , with increases in average whole body weight of over 38% compared to unselected bay scallops, have been selected from the F 1 hybrids between the Peruvian scallop and A. i. irradians 6,7 . Recently, a new strain was www.nature.com/scientificdata www.nature.com/scientificdata/ obtained by further crossing the 'Bohai Red' strain with A. i. concentricus. Interestingly, this strain exhibited a better tolerance to high temperature than the 'Bohai Red' strain and had a longer life span than A. i. concentricus (Zhigang Liu, personal communication). In addition, the selection of a genetically stable strain in bivalves by traditional breeding methods could take six to ten years, but marker-assisted selection based on genomic data can greatly reduce breeding duration.
In addition to its application in breeding, genomic data can also be immediately employed in studies of evolution, adaptation, longevity, gonad development, and sex determination in bivalves [8][9][10][11] . To date, several genomes have been sequenced and assembled in bivalves. For example, assembly of the Pacific oyster (Crassostrea gigas) genome provided insights into how sessile oysters adapt to adverse environments 8 . Analyses of the Japanese scallop (Patinopecten yessoensis) genome, the first sequenced scallop genome, revealed that scallops may have a conserved primitive karyotype close to that of the ancestral bilaterian 9 . Comparison of the genome sequences of a deep-sea mussel (Bathymodiolus platifrons) and a shallow-water mussel (Modiolus philippinarum) deepened our understanding of how deep-sea organisms adapt to extreme environments 10 .
In a previous study, we have sequenced and assembled the genome of the Peruvian scallop 12 . The genomic data of Argopecten scallops and their hybrids will allow us to investigate the evolutionary relationships among Argopecten scallop species and subspecies, to study the molecular mechanisms underlying scallop adaptations to diverse habitats and to understand their wide variation in life span as well as the development of male sterile gonads in their F 1 hybrids.
In this study, we sequenced and assembled the genomes of the two bay scallop subspecies, A. i. irradians and A. i. concentricus. Together with the genomic data of the Peruvian scallop, the results of our present study and subsequent genome-wide association studies will eventually facilitate the breeding progress in these Argopecten scallops.

Methods
Sample collection and genomic DNa extraction. Genomic DNA was extracted from the adductor muscle of a single specimen from a pure line of A. i. irradians (Fig. 1a) and A. i. concentricus (Fig. 1b), which have been bred by self-fertilization in a scallop farm in Laizhou, Shandong Province, China. The quality of the DNA samples was checked by electrophoresis on 1% agarose gels. The purity of the DNA was also checked using a NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA). DNA concentration was measured using a Qubit ® DNA Assay Kit in Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA). A total of 1.5 μg DNA per sample was used for subsequent sample preparations.
DNa library preparation and whole genome sequencing. Sequencing libraries were generated using the Truseq Nano DNA HT Sample Preparation Kit (Illumina, USA) following the manufacturer's instructions. Index codes were used to cross-index the sequences and samples, that is, the DNA samples were fragmented by sonication and then end-polished, A-tailed, and ligated with the full-length adapter for Illumina sequencing followed by PCR amplification. The resulting PCR products were purified (AMPure XP system) and the sequence libraries were analyzed for size distribution by Agilent2100 Bioanalyzer and quantified using real-time PCR.
These libraries were sequenced on an Illumina HiSeq4000 platform using a 150-bp paired-end sequencing protocol. Additional sequencing was performed on a PacBio Sequel instrument with a Sequel TM Sequencing Kit 1.2.1 (Pacific Biosciences, USA) for both subspecies. Raw reads in the FASTQ format were first processed using Trimmomatic. In this step, clean reads were obtained by removing reads containing adapter sequences, poly-N repeats, and low-quality sequences. In addition, the Q20, Q30, and GC content of the clean reads was determined for quality control. All downstream analyses were based on the high quality, clean reads identified here.
Genome assembly. To estimate the genome size of both subspecies, a routine 17-mer frequency distribution analysis 13 was performed according to the following formula: genome size = k-mer number/peak depth (Table 1). A routine assembly strategy was applied for the genome assemblies of both scallops. Briefly, all high-quality reads were assembled into scaffolds using Platanus v1.2.4 14 , and the gaps were filled using GapCloser 15 . Subsequently, PBJelly v14.1 was applied for additional gap filling with Pacbio reads 16 . Finally, all Illumina reads were employed to correct the genome assemblies in Pilon v1.18 for two rounds 17 (Table 1). Genome assessment. Following the initial assembly, the integrity of both genome assemblies was assessed by mapping the reads from short-insert libraries onto the assembled genomes using Burrows-Wheeler Aligner www.nature.com/scientificdata www.nature.com/scientificdata/ (BWA) 18 , which can align the clean reads from multiple samples against the reference genomes (settings: bwa mem -t 4 -k 32 -M -R). Alignment files were converted to BAM files using SAMtools 19 . In addition, potential PCR duplications were removed using the SAMtools with command "rmdup". If multiple read pairs had identical external coordinates, only the pair with the highest mapping quality was retained. Subsequently, the Core Eukaryotic Genes Mapping Approach (CEGMA) was employed to evaluate the completeness of both genome assemblies 20 . Among the 248 ultraconserved core eukaryotic genes (CEGs), we identified 231 (93.15%, complete + partial) and 227 (91.53%, complete + partial) CEGs in the genomes of northern and southern subspecies, respectively. Benchmarking Universal Single-Copy Orthologues (BUSCO) was used to quantitatively assess the completeness of genome assembly based on evolutionarily-informed expectations of gene content from near-universal single-copy orthologues 21 . The assessment demonstrated that 91% of the 843 single-copy genes were identified in both genome assemblies, containing C: repeat annotation. Two methods were employed to identify transposable elements (TEs) in the assembled genomes. When using the ab-initio method, RepeatModeler was used to build a species-specific repeat database (parameters set as '-engine_db wublast') 22 . When using a homology-based method, RepeatMasker 23 was employed to identify repeats with known repeat libraries (Repbase) 24 using the following parameters: '-a -nolow -no_is -norna -parallel 3 -e wublast-pvalue 0.0001' , along with RepeatProteinMask (the parameter set as '-noLowSimple -pvalue 0.0001 -engine wublast') 23 , and the repbase data were collected from a comprehensive database of undifferentiated species (RepBase Metadata Database RELEASE 20170127). In addition, tandem repeats were identified using Tandem Repeats Finder (TRF) with the parameters setting as 'Match = 2, Mismatching penalty = 7, Delta = 7, PM = 80, PI = 10, Minscore = 50, MaxPeriod = 2,000' 25 (Table 2). Gene annotation. de novo prediction. Protein-coding genes in the assembled genomes were annotated using de novo prediction by homology with transcriptome data-based evidence. Four programs were employed for the de novo prediction of genes, including Augustus v3.2.1 (with the following parameters: '-uniqueGeneId true -noInFrameStop = true -gff3 on -genemodel complete -strand both') 26 , Genscan (using default parameters) 27 , GlimmerHMM (with the following parameters: ' -f -g') 28 , and SNAP (using default parameters) 29 .  www.nature.com/scientificdata www.nature.com/scientificdata/ purpuratus) were used for homology-based searches against the two genome assemblies using TBLASTn (e-value ≤ 10 −5 ) 30 . The final gene structures were predicted using GeneWise (with the following parameters: '-genesf ') 31 .
Transcriptome-based annotation. Transcriptome data from different tissues including kidney, hepatopancreas, and haemolymph were mapped onto each genome assembly using Tophat (with the following parameters: 'max-intron-length 500000 -m 2-library-type fr-unstranded') 32 , and used for gene modeling using Cufflinks (with the following parameters: '-multi-read-correct') 33 according to the pair-end relationships and the overlaps between aligned reads.
Gene set integration. Following de novo prediction, homology-based annotation, and transcriptome-based prediction, we integrated the gene models using EvidenceModeler (EVM) 30 to generate a comprehensive and non-redundant gene set (Table 1).
Functional assignment. Gene function annotation was performed by aligning the predicted protein sequences against various protein databases-including the SwissProt 34 and NCBI non-redundant (Nr) databasesusing BLASTP (e-value ≤ 10 −5 ). Gene domain annotation was performed by searching the InterPro database 35 . All genes were aligned against the Kyoto Encyclopedia of Genes and Genomes (KEGG) 36 database to identify gene pathways. Gene Ontology (GO) terms of the genes were obtained from the corresponding InterPro entry 37 . ortholog and gene family expansion analysis. The protein-coding genes from both scallop genome assemblies and seventeen other sequenced species including Pacific oyster, owl limpet, Amphioxus, nematode, fruit fly, leech, human, octopus, red flour beetle (Tribolium castaneum), polychaete (Capitella teleta), brachiopod (Lingula anatina), sea slug (Aplysia californica), abalone (Haliotis discus), pearl oyster (Pinctada fucata), Yesso scallop (P. yessoensis), deep-sea vent mussel (B. platifrons) and shallow-water mussel (M. philippinarum) were analyzed. All data were downloaded from the Ensembl 38 or NCBI 39 databases. Gene family analysis was performed based on the homologs of the protein-coding genes in the related species, which was initially implemented by the alignment of an "all against all" BLASTP. Subsequently, alignments with high-scoring segment pairs (HSPs) were conjoined for each gene pair by Solar 40 to process the mapping results. To identify homologous gene-pairs, more than 30% coverage of the aligned regions in both homologous genes was required. Finally, homologous genes were clustered into gene families by OrthoMCL 41 . A p-value cut-off of 1e-5 was chosen for putative orthologues or paralogs, which were converted into a graph for the nodes of representative protein sequences. The resulting graph is represented by a symmetric similarity matrix to which an MCL algorithm was applied (with the following parameters: "-inflation 1.5") to regulate cluster tightness (Fig. 2a,b).

Genome evolution analysis. Phylogenetic analysis. To trace the evolutionary position of A. i. irradians
and A. i. concentricus, a dataset containing 107 single-copy protein-coding genes retrieved from the 19 species mentioned above was used for phylogenetic tree construction and divergence time estimation. Protein sequences for these single-copy genes were aligned by MUSCLE 42 one by one, and then were concatenated to the final dataset. ProtTest 43 was used to select the best-fit model for this dataset. Then, the phylogenetic tree was reconstructed using the RAxML method (version 7.2.3) 44 with LG + G + I + F model with the proportion of invariable sites 0.07 and Gamma shape parameter 0.83. The clade containing H. sapiens and B. floridae was set as outgroup. Clade support was assessed using the bootstrapping algorithm in the RAxML with 1000 alignment replicates.
Estimation of divergence time. Species divergence time was inferred based on the same dataset containing 107 single-copy protein-coding genes from the 19 species using the MCMCTree function included in PAML v4.7a 45 with the following parameters: '-model 0-rootage 1200 -clock 3' . For their divergence time estimation, reference divergence times obtained from TimeTree database 46

Data records
The whole genome sequences of A. i. irradians and A. i. concentricus were deposited in public repositories. The raw sequencing and transcriptomic data were deposited in NCBI Sequence Read Archive, under the SRA study accession SRP174526 47 . This whole-genome project including the assembly fasta, annotation and protein sequencing was uploaded to Dryad (https://doi.org/10.5061/dryad.hdr7sqvdr) 48 . All genome annotation and phylogenetic tree files were uploaded to Figshare (https://doi.org/10.6084/m9.figshare.c.4356239) 49 . The genome assemblies are also available at the NCBI Assembly website 50,51 .

technical Validation
To produce high-quality draft genome assemblies, we applied whole-genome sequencing, assembly, and annotation of the two bay scallop subspecies. The whole genome shotgun sequencing strategy was used for both bay scallop subspecies. We constructed six sequencing libraries including two short-insert libraries (250 bp and 450 bp) and four long-insert libraries (2, 5, 10, and 20 kb) for A. i. irradians and A. i. concentricus, respectively. For A. i. www.nature.com/scientificdata www.nature.com/scientificdata/ irradians, a total of 253.17 gigabases (Gb) of raw reads were generated while a total of 272.97 Gb of raw reads were generated for A. i. concentricus. For A. i. irradians, a total of 3.86 × 10 10 k-mers with a peak k-mer depth of 38 were employed to obtain an estimated genome size of 996.07 Mb (Table 1). In A. i. concentricus, a total of 4.97 × 10 10 k-mers and a peak k-mer depth of 50 were employed to obtain the estimated genome size of 974. 3 Fig. 3 Phylogenetic position of the sequenced species. The phylogenetic tree was constructed based on a dataset from 107 single-copy orthologues using the RAxML method. Clade support was assessed using the bootstrapping algorithm with 1,000 alignment replicates. (a) The phylogenetic tree was reconstructed using the RAxML method with LG + G + I + F model. The tree is drawn to scale, with branch lengths proportional to the number of amino acid substitutions. Bootstrap values are presented above the nodes. (b) Species divergence time was estimated using the MCMCTree function in the PAML with the parameter of '-model 0-rootage 1200 -clock 3' . Red nodes in the phylogenetic tree represented the reference divergence times, which were applied to calibrate the divergence dates of these examined species.
www.nature.com/scientificdata www.nature.com/scientificdata/ with a coverage of 86.41%. These mapping results suggest good reliability for both genome assemblies, which are close to the assembly of the Peruvian scallop genome in our previous study but better than those of other related bivalve species 12 .
A protein is classified as complete if the alignment of the predicted protein to the HMM profile represents at least 70% of the original KOG domain, otherwise, it is classified as partial. Our evaluation results demonstrated that both genome assemblies covered 231 (93.15%) of the 248 Core Eukaryotic Gene sequences, indicating a high level of completeness within the two genome assemblies. A related assessment identified 91% of the 843 single-copy genes in both genome assemblies. These data indicate the high integrity of both genome assemblies. A total of 389,681,429 and 412,788,948 bp of repeat sequences were predicted in the A. i. irradians and the A. i. concentricus genomes, respectively. These repeat sequences accounted for 46.43% and 47.17% of the corresponding genome assemblies (Table 2). A total of 26,777 protein-coding genes were predicted in A. i. irradians with an average transcript length of 11.86 kb. The public functional databases Swissport, interpro and NR were used for gene prediction and annotation. Similarly, a total of 25,979 protein-coding genes were predicted in A. i. concentricus with an average transcript length of 12.17 kb (Table 1). In total, 24,943 (93.2%) and 24,428 (94%) predicted proteins could be functionally annotated in A. i. irradians and A. i. concentricus, respectively, using public databases. In total, the protein-coding genes were classified into 48,052 gene families and 107 strict single-copy orthologs (Fig. 2a). Compared to other examined species, 1,949 and 1,769 gene families were exclusively presented in A. i. irradians and A. i. concentricus, respectively (Fig. 2b).
The evolutionary position and divergence time of A. i. irradians and A. i. concentricus were elucidated in this study. The results of the phylogenetic tree showed that outgroup clade containing H. sapiens and B. floridae located in the basal position of the whole tree with high confidence (bootstrap value = 100%). Meanwhile, we found that A. i. irradians and A. i. concentricus clustered together with the 100% bootstrap value (Fig. 3a), and then merged as a sister group to P. yessoensis, as it did in the divergence time tree (Fig. 3b). It shows that this clade has a close relationship with the other two clades containing C. gigas, P. fucata and M. philippinarum, B. platifrons. Besides, we estimated the divergence times of A. i. irradians and A. i. concentricus using single-copy protein-coding genes from the 19 examined species (Fig. 3b). The result showed that the divergence time between the Northern subspecies (A. i. irradians) and the Southern subspecies (A. i. concentricus) happened at ~26.4 Mya ago, and the analysis suggested that the ancestor of A. i. irradians, A. i. concentricus and P. yessoensis originated ~85.9 Mya.

Code availability
In the study, we did not use any custom specific code. The command line for each step is executed as indicated for each step of all bioinformatics procedures.