Sequencing and quality evaluation

A total of 639.75 Gb raw reads were generated for S. leurokolos. After filtering and correction, a total of 599.63 Gb clean reads were derived (Table 1). The Q20 and Q30 values of the sequencing data were 96.28 and 91.18%, respectively (Table 1). It has been specified that Q20 and Q30 values should be at least 90 and 85% (Li et al., Reference Li, Song, Jin, Li, Gong and Wang2019). Therefore, the sequencing data of S. leurokolos genome show extreme precision in the present study. GC content is an important factor in many experiments and bioinformatic analysis, especially for next-generation sequencing where the sequenced DNA has gone through multiple rounds of PCR amplification. High or low GC content will reduce sequencing coverage and cause sequencing bias (Bentley et al., Reference Bentley, Balasubramanian, Swerdlow, Smith, Milton, Brown, Hall, Evers, Barnes and Bignell2008; Aird et al., Reference Aird, Ross, Chen, Danielsson, Fennell, Russ, Jaffe, Nusbaum and Gnirke2011; Cheung et al., Reference Cheung, Down, Latorre and Ahringer2011). In this study, GC content of S. leurokolos sequences was 37.6% showing a mid GC content (30–47%) (Shangguan et al., Reference Shangguan, Han, Kayesh, Sun, Zhang, Pervaiz, Wen and Fang2013). Overall, these results indicate high-quality sequencing data obtained for S. leurokolos.

Table 1. Summary information for the S. leurokolos genome sequencing and genome assembly

Q20: the ratio of data with accuracy above 99% in total data. Q30: the ratio of data with accuracy above 99.90% in total data

K-mer analysis and genome size estimation

The genome size, heterozygosity and repetitive ratio of S. leurokolos were evaluated using K-mer distribution analysis, and the 17-mer yielded the highest model fit (Figure 2 and Table 2). K-mer analysis revealed a unique bimodal profile with a high heterozygous peak around 50× coverage and a lower homozygous peak around 100× coverage (Figure 2). By calculation, the genome size of S. leurokolos was estimated to be 5.08 Gb (Table 2). Flow cytometry is another method for the prediction of genome size. Previous study for other four alvinocaridid shrimps based on flow cytometry reveals that genome sizes range from 10,160 Mp in Rimicaris exoculata to 13,050 Mp in Chorocaris chacei (Bonnivard et al., Reference Bonnivard, Catrice, Ravaux, Brown and Higuet2009), displaying a large genome size in the family Alvinocarididae. It seems that the genome size of S. leurokolos is much smaller than those of other alvinocaridid shrimps, or its genome size has been underestimated by GSS. The significant discordance of genome size revealed by GSS and flow cytometry has been also detected in other decapods, such as crayfish Procambarus clarkia, showing larger genome size by flow cytometry analysis than that revealed by GSS (Shi et al., Reference Shi, Yi and Li2018). However, muscle instead of haemolymph cell has been used in the flow cytometry analysis for alvinocaridid shrimps (Bonnivard et al., Reference Bonnivard, Catrice, Ravaux, Brown and Higuet2009), which might be due to the difficulty in collecting living shrimp samples from deep sea. It may influence the quality of cell suspension preparation and in turn affect the precision of genome size estimation. On the other hand, the high heterozygosity and repeat ratio characteristics of S. leurokolos genome as shown below might bring biased results in genome size estimation by affecting the K-mer depth distribution (Shi et al., Reference Shi, Yi and Li2018). In brief, GSS and flow cytometry should be combined to estimate genome sizes of deep-sea species with large and complex genome, and the genome size of S. leurokolos might be larger than 5.08 Gb.

Figure 2. K-mer (K = 17) analysis for estimation of the genome size of S. leurokolos. The x-axis represents coverage, and the y-axis represents the frequency at each depth.

Table 2. Statistics of the estimated S. leurokolos genome size and other characteristics

According to the K-mer distribution, an extremely high heterozygosity 2.85% was detected in S. leurokolos genome (Figure 2 and Table 2). It has been suggested that genome assembly will be difficult if the heterozygosity rate exceeds 0.5%, and it is even more difficult if it exceeds 1% (Marçais and Kingsford, Reference Marçais and Kingsford2011). The repeat ratio of S. leurokolos genomic sequences was also high (87.03%) (Figure 2 and Table 2). The high heterozygosity rate and repeat ratio have been also revealed in other decapods, such as Litopenaeus vannamei, Penaeus chinensis and P. monodon (Zhang et al., Reference Zhang, Yuan, Sun, Li, Gao, Yu, Liu, Wang, Lv, Zhang, Ma, Wang, Lin, Wang, Zhu, Zhang, Zhang, Jin, Yu, Kong, Xu, Chen, Zhang, Sorgeloos, Sagi, Alcivar-Warren, Liu, Wang, Ruan, Chu, Liu, Li and Xiang2019; Van Quyen et al., Reference Van Quyen, Gan, Lee, Nguyen, Nguyen, Tran, Nguyen, Khang and Austin2020; Uengwetwanit et al., Reference Uengwetwanit, Pootakham, Nookaew, Sonthirod, Angthong, Sittikankaew, Rungrassamee, Arayamethakorn, Wongsurawat, Jenjaroenpun, Sangsrakru, Leelatanawit, Khudet, Koehorst, Schaap, Martins dos Santos, Tangy and Karoonuthaisiri2021; Yuan et al., Reference Yuan, Zhang, Wang, Sun, Liu, Li, Yu, Gao, Liu, Zhang, Kong, Fan, Zhang, Feng, Xiang and Li2021b; Wang et al., Reference Wang, Ren, Liu, Li, Lv, Wang, Zhang, Wei, Zhou, He and Li2022b), and difficulties in genome assembly seem to be common problem in decapods due to high heterozygosity and repeat ratio (Yuan et al., Reference Yuan, Zhang, Li and Xiang2021a).

Genome de novo assembly

To assemble the draft genome of S. leurokolos, two K-mer values, 41 and 63 bp were selected. Unfortunately, too much computer memory was required and the assembly task could not be completed when using the 41 bp K-mer value. A complete assembly using 63 bp K-mer value was obtained (Table 1). Finally, our efforts recovered a total of 9,527,856,577 bp scaffolds with the scaffold N50 value of 597 bp, and the maximum scaffold was 69,344 bp in length (Table 1). It is apparent that the size of draft genome assembly is almost twice as large as the estimated genome size based on 17-mer analysis. The most plausible explanation for the genome assembly size deviation may be that the presence of a large number of repetitive elements (87.03%) and high heterozygosity (2.85%) of S. leurokolo genome might induce the assembly has multiple copies of the same genomic region and even contained two divergent haplotypes (Pflug et al., Reference Pflug, Holmes, Burrus, Johnston and Maddison2020; Hu et al., Reference Hu, Feng, Xiang, Wang, Salojärvi, Liu, Wu, Zhang, Liang, Jiang, Liu, Ou, Li, Fan, Mai, Chen, Zhang, Zheng, Zhang, Peng, Yao, Wai, Luo, Fu J, Tang, Lan, Lai, Sun, Wei, Li, Chen, Huang, Yan, Liu, McHale, Rolling, Guyot, Sankoff, Zheng, Albert, Ming, Chen, Xia and Li2022; Wyngaard et al., Reference Wyngaard, Skern-Mauritzen, Malde, Prendergast and Peruzzi2022). The average GC content of S. leurokolos assembled genome was about 36.12%. To further evaluate the data of our assembly, we compared it to previously reported genome survey data of decapods. The scaffold N50 of S. leurokolos is much shorter than that of Pacific white shrimp L. vannamei (1343 bp) (Yu et al., Reference Yu, Zhang, Yuan, Li, Chen, Zhao, Huang, Zheng and Xiang2015) and red swamp crayfish P. clarkia (1426 bp) (Shi et al., Reference Shi, Yi and Li2018). The inherent defects of second-generation sequencing technology in read length and high complexity of the large genome of S. leurokolos itself should be the main reasons for the poor assembly. We hold the opinion that the large and complex genome of S. leurokolos represents typical challenges faced by all alvinocaridid shrimp genomes, which partly explains why genomic resources for alvinocaridid shrimps are so limited compared to those of many other deep-sea organisms. Hence, developing new assemblers and bioinformatics tools and using combination of short- and long-read sequencing technologies (i.e. PacBio, Oxford Nanopore Technologies, ONT) are expected to solve these challenges for assembling a high-quality genome. The current GSS data could serve as a reference for subsequent whole-genome sequencing project of S. leurokolos.

Genomic repetitive elements annotation

Repetitive sequences, especially transposable elements (TEs), are known to be an evolutionary precursor of many genes, a driving force in the evolution of epigenetic regulation and an important factor in genomic stability maintenance and evolution (Jurka et al., Reference Jurka, Kapitonov, Kohany and Jurka2007). In total, 4250 Mb repetitive elements were identified in S. leurokolos draft genome, accounting for 44.62% of the assembled genome (Table 3). Combining the results from RepeatMasker and RepeatProteinMask analyses, our results revealed that among these repetitive sequences, 38.92% (3708 Mb) were TEs, but 16.49% could not be classified within the TEs (Table 4). Long interspersed nuclear elements (LINEs) were the most common among the TEs, accounting for 10.45%, followed by DNA transposons (6.09%) and long-terminal repeat elements (LTRs) (4.79%) (Table 4). These repetitive elements, including LINEs, DNA and LTRs, also take up a large proportion of genomes in many other decapod crustaceans (Baeza, Reference Baeza2020; Tang et al., Reference Tang, Wang, Liu, Zhang, Jiang, Li, Wang, Sun, Sha, Jiang, Wu, Ren, Li, Xuan, Ge, Jiang, She, Sun, Qiu, Wang, Wang, Qiu, Zhang and Li2020; Chak et al., Reference Chak, Harris, Hultgren, Jeffery and Rubenstein2021; Uengwetwanit et al., Reference Uengwetwanit, Pootakham, Nookaew, Sonthirod, Angthong, Sittikankaew, Rungrassamee, Arayamethakorn, Wongsurawat, Jenjaroenpun, Sangsrakru, Leelatanawit, Khudet, Koehorst, Schaap, Martins dos Santos, Tangy and Karoonuthaisiri2021). However, it has been suggested that the ‘unclassified’ TEs with a large proportion may contain species-specific variants of known repetitive elements, and we should be cautious when comparing these datasets directly with those of other species (Murgarella et al., Reference Murgarella, Puiu, Novoa, Figueras, Posada and Canchaya2016).

Table 3. Statistics of repetitive sequence annotation in the S. leurokolos draft genome assembly

Table 4. Statistics of TEs in the S. leurokolos draft genome assembly

RepBase TEs and TE proteins were obtained based on the RepBase library using RepeatMasker and RepeatProteinMask, respectively. De novo repeat prediction was performed using RepeatMasker against the de novo repeat library of S. leurokolos, which was constructed by the programs LTR_FINDER and RepeatModeler. Combined TEs were the union of the three methods.

Microsatellite analysis

It is widely recognized that as a most popular and versatile genetic marker, SSRs are widely used for the genetic characterization of populations due to their abundance in genome, high polymorphism and co-dominant nature (Abdul-Muneer, Reference Abdul-Muneer2014). In the assembled scaffolds, a total of 12,121,553 microsatellite motifs were identified in S. leurokolos (Table 5). Among them, the di-nucleotide was the most abundant, accounting for 70.27% of the total SSRs, which was followed by tri- (25.54%), tetra- (3.33%), penta- (0.50%) and hexa- (3.36%) nucleotide SSRs (Table 6). Our finding shows that both di-nucleotide and tri-nucleotide SSRs are numerous, and the number of repetitions is inversely proportional to the length of repetitions. This result is consistent with those in other crustaceans, such as kuruma prawn Marsupenaeus japonicus (Lu et al., Reference Lu, Luan, Kong, Hu, Mao and Zhong2017), Japanese mantis shrimp Oratosquilla oratoria (Cheng et al., Reference Cheng, Zhang and Sha2018) and Antarctic krill Euphausia superba (Huang et al., Reference Huang, Bian, Liu, Wang, Xue, Huang, Yi, You, Song, Mao, Song and Shi2020). It has been proposed that longer repeats have downward mutation bias and short persistence times (Harr and Schlötterer, Reference Harr and Schlötterer2000), and therefore, less SSRs with longer repeat units exist in genomes.

Table 5. Statistics of SSR distribution in the S. leurokolos draft genome assembly

Table 6. Organization of the S. leurokolos mitogenome

Mitochondrial genome and candidate molecular marker identification

Mitochondria are essential organelles that generate most chemical energy to power the cell's biochemical reactions. There is evidence that mitochondrial DNA plays a role in many aspects of biological life history, such as lifespan, fertility, resistance to starvation, altitude adaptation and regulation of temperature (Ballard and Melvin, Reference Ballard and Melvin2010). It is therefore of significant importance to investigate the mitochondrial genome of S. leurokolos inhabiting deep-sea chemosynthetic ecosystems. In this study, we assembled a 15,906 bp long complete mitochondrial genome (GenBank accession no. OQ622002) of S. leurokolos from the GSS data. It consisted of 13 protein-coding genes (PCGs), 2 ribosomal RNA genes (rrnS and rrnL), 22 transfer (tRNA) genes and a non-coding hypervariable control region (1026 bp) between rrnS and tRNA-Ile, showing the typical alvinocaridid shrimp mitogenome arrangement model (Table 6). Most of the PCGs and tRNA genes were encoded on the positive strand. Gene overlaps in 19 gene junctions (a total of 57 bp in length) and intergenic spaces in 14 gene junctions (ranging from 1 to 50 bp) were also observed (Table 6).

Moreover, mitochondrial DNA fragments have been proved to be efficient molecular markers in phylogenetic and population genetic analysis. In order to identify candidate markers, we aligned the mitochondrial genome assembled in this study with the previous reported S. leurokolos mitochondrial genome (Sun et al., Reference Sun, Hui, Wang and Sha2018a). By comparison, 3 indels (all located in the control region) and 71 SNPs were detected. The SNPs included 66 transitions and 5 transversions: 47 in PCGs, 3 in tRNAs, 1 in rRNAs and 19 in non-coding regions. Of the 47 SNPs in PCGs, only four mutations were non-synonymous substitutions (Table 7), which occurred in cox1, nad2, cytb and nad1 (Table 7). It is a general observation in molecular evolution that functional importance and substitution rate are negatively correlated (Sun et al., Reference Sun, Li and Kong2010). This means that the more functionally important genes (or genetic regions) evolve more slowly due to their important effects or strong functional constraints (Kimura, Reference Kimura1983; Yang, Reference Yang2006). In addition, the relatively high substitution rates observed in tRNA-Ala (1.59%), control region (1.58%), tRNA-Cys (1.49%) and tRNA-Trp (1.33%) may indicate relatively low functional constraints in these regions.

Table 7. Summary of SNPs in S. leurokolos mitochondrial genome

To date, population genetic and phylogenetic studies for alvinocaridid shrimps are mainly based on mitochondrial cox1, 12S rDNA and 16S rDNA genes (Yahagi et al., Reference Yahagi, Watanabe, Ishibashi and Kojima2015; Sun et al., Reference Sun, Sha and Wang2018b). In this study, cox1, nad2, nad4 and control region show high mutation rate, and the sequences are long enough for primer design. Hence, these mitochondrial genes can be selected as candidate markers for population genetic studies for S. leurokolos. However, it requires further validation by amplification and sequencing in more individuals.