A chromosome-level genome assembly of the Asian giant softshell turtle Pelochelys cantorii

The Asian giant softshell turtle Pelochelys cantorii is one of the largest aquatic turtles in China and has been designated a First Grade Protected Animal in China. To advance conservation research, a combination of Illumina short-read, PacBio long-read, and Hi-C scaffolding technologies was used to develop a high-quality chromosome-level genome assembly for P. cantorii. A total of 262.77 Gb of clean data were produced (121.6 × depth) and then the genome was assembled into 2.16 Gb with a contig N50 of 41.44 Mb and scaffold N50 length of 120.17 Mb, respectively. Moreover, about 99.98% assembly genome sequences were clustered and ordered onto 33 pseudochromosomes. Genome annotation revealed that 21,833 protein-coding genes were predicted, and 96.40% of them were annotated. This new chromosome-level assembly will be an enabling resource for genetic and genomic studies to support fundamental insight into P. cantorii biology.


Background & Summary
The Asian giant softshell turtle Pelochelys cantorii is one of the largest aquatic turtles and is widely distributed in Southeast Asia.With rapid economic development, the numbers of P. cantorii have sharply declined due to overhunting and habitat destruction in China.Although listed as a national first-class protected animal in China in 1989, wild P. cantorii individuals in China have been reported at a declining rate in the past 30 years 1,2 , and only 13 P. cantorii individuals from the wild are in captivity across 6 different locations.P. cantorii is extremely endangered in China 1,3 .The fate of another large soft-shelled turtle in China, the Yangtze giant softshell turtle, Rafetus swinhoei, is even more worrisome.No individuals were found in the wild, and two individuals (1 female and 1 male) reared in captivity failed in assisted reproduction after 2008.In 2019, the lone female R. swinhoei died during artificial insemination.If there are no R. swinhoei individuals in the wild, the fate of R. swinhoei is sealed.If the protection of P. cantorii is not strengthened, it is inevitable that the P. cantorii will not become the next R. swinhoei.
In 2014, 10 cantorii hatchlings were successfully bred in captivity from 2 turtle individuals (1 female and 1 male) 4 .From 2015 to 2021, four captive P. cantorii individuals (2 females and 2 males) successfully bred and raised 950 juveniles that are currently 1-7 years old.Knowledge of the artificial breeding of P. cantorii offers hope for the conservation of this species 5 .With the breakthrough of the artificial breeding of P. cantorii, related protection work has received attention from the government of China.P. cantorii is listed as one of the nine aquatic wild animals under the key management of the Ministry of Agriculture and Rural Affairs (MARA) of the People's Republic of China.In 2019, the MARA of China issued the "Pelochelys cantorii Rescue Action Plan (2019-2035)".In September 2020, the MARA of China organized and carried out the first wild adaptation protection test for P. cantorii.A total of 20 juvenile turtles with implanted PIT chips were released in a reservoir in Gaoming, Foshan.These turtles were between 4-5 years old and weigh 1.04-1.66kg.
With the development of sequencing technology, conservation genomics was born, which overcomes the limitations of conservation genetics, such as the lack of markers to a large extent, and helps to solve some unresolved problems in the conservation biology of organisms, especially endangered species 6 .For example, de Novo sequencing and population resequencing analysis based on Yangtze finless porpoise Neophocaena Asiaeorientalis asiaeorientalis more accurately determined the phylogenetic relationships and population genetic structure, and confirmed the independent species status of the Yangtze finless porpoise 7 .Through the analysis of genomic data of hot spring snakes, the genes involved in DNA damage repair (FEN1) and hypoxia response (EPAS1) were identified, providing research ideas for the mechanism of environmental adaptation.With the continuous development and optimization of high-throughput sequencing technology, research on conservation genomes will continue to improve.
Third-generation sequencing technologies, such as Pacific Biosciences (PacBio) Single Molecule Real-Time (SMRT) and Oxford Nanopore sequencing, characterized by single-molecule sequencing, can overcome the shortcomings of second-generation sequencing methods 8 .Hi-C, based on the combination of chromatin conformation capture and high-throughput sequencing, reveals the interaction information between chromosome fragments through the analysis of sequencing data [9][10][11][12] .
Tortoises are an ancient reptile with more than 300 species and have significant evolutionary and ecological value.No more than 10 turtle genome projects have been completed since the first investigation on turtle genomes were published in 2013 (Chinese soft-shell turtle and green turtle).Therefore, it is urgent to expand the diversity of the turtle genome databases.
Here, we sequenced and de novo assembled the genome of a young P. cantorii via a combination of long-read PacBio Sequel II platform and Hi-C sequencing technology.The high-quality reference genome constructed in this study will not only be of benefit to serve as the genetic basis for in-depth investigations of turtle evolution and biology but also offers a valuable genetic resource for turtle conservation.

Methods ethics statement.
All the experimental procedures regarding the turtle involved in this experiment were approved by the Experimental Animal Care and Ethics Committee of the Pearl River Fisheries Research Institute, Chinese Academy of Fishery Sciences.

Sample preparation and genome sequencing.
A healthy 2-year-old adult was obtained from the breeding center of the Pearl River Fisheries Research Institute, Guangzhou, Guangdong, China (Fig. 1).The turtle was anaesthetized using tricaine methanesulfonate (MS-222, Sigma) before sampling.The fresh muscle tissue was cut to small sizes (about 20 µg) and conducted for DNA extraction using a DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA, USA), according to the manufacturer's instructions.Only high-quality DNA can be used for the following library preparation, sequencing, and Hi-C library construction, so the quantity and quality of the genomic DNA were first evaluated by using Qubit 3.0 (Life Invitrogen, USA) and 1% agarose gel electrophoresis, respectively.The library was then constructed and sequenced on Illumina HiSeq4000 (Illumina, San Diego, CA, USA) and the PacBio Sequel II platform according to the manufacturer's instructions.
Tissues from the heart, liver, kidney, muscle, eye, blood, and lung of the same turtle were used for RNA extraction by TRIzol Universal Reagent (TIANGEN Biotech, Beijing, China).Two micrograms of total RNA from each tissue were pooled for RNA sequencing on a PacBio Sequel II.
Genome size estimation and initial genome assembly.Genome size was estimated by k-mer analysis 13 using filtered reads from the three 350 bp libraries constructed using genomic DNA and a k-mer size of 21.Genome size, repetition ratio, and heterozygosity were assessed by k-mer depth frequency distribution analysis (Table 1).The average k-mer depth was 56, and the total number of k-mers obtained from sequencing data was 131,112,165,086, further estimating the genome size of P. cantorii (Fig. 2 and Table 1).Sequences with a k-mer depth greater than 113 were repeated sequences, and those with a k-mer depth of approximately 28 were heterozygous sequences.After removing k-mers of abnormal depth, the genome size of P. cantorii was evaluated to be 2.20 Gb (based on the following formula: G = k-mer number/mean k-mer depth), with a low heterozygosity of approximately 0.14% and a repeat sequence content of 26.29% (Table 1).The PacBio sequencing platform generated a total of 262.77Gb of clean data (approximately 121.60 × ), with a read N50 of 22.85 kb and an average read length of 14.62 kb (Table 2).The data were first assembled with WTDBG2, which generated a genome assembly with a total length of 2.16 Gb, 548 contigs, and 41.44 Mb contig N50 (Table 3).

Chromosome assembly using Hi-C data.
A high-throughput chromatin conformation capture (Hi-C) library was then created for sequencing from 2-year-old individual muscle tissue.The sample was fixed with paraformaldehyde and enzymatically digested with DnpII, generating sticky ends.The sticky ends can be repaired by ' A' or 'C' deoxynucleotides with biotin under the action of DNA polymerase.After labeling the 5′ overhang with a biotinylated residue, the DNA fragments were ligated to each other to form chimeric circles using DNA ligase.Then, the ligated DNAs were uncrosslinked, purified, and sheared to 300 bp-700 bp fragments, followed by capture with streptavidin beads and preparation for sequencing.Finally, the Hi-C libraries were quantified and sequenced using the Illumina HISeq X 10 platform in PE150 mode with paired-end methods.To obtain high-quality Hi-C data, the joint information, low-quality bases, and undetected bases in the original sequencing data were first removed to generate clean reads.Second, BWA v0.7.10-r789 was applied to align the clean reads to the draft genome assembly with the default parameters.HiC-Pro v2.invalid read pairs, including dangling-end and self-cycle, re-ligation, and dumped products.The contigs of the primary genome assembly were corrected by splitting them into segments of 50 kb on average.After checking any two segments that showed inconsistent connection with information from the raw scaffold, LACHESIS (ligating adjacent chromatin enables scaffolding in situ) was used to produce chromosome-level scaffolds.The final assembled genome was 2.16 Gb and anchored 532 contigs.The scaffold N50 was 41.44 Mb with a fixed rate of 95.86% (Table 3).Among the sequences located on chromosomes, the length of the sequence that can determine the order and direction is up to 99.97%.The Hi-C scaffolding of the genome resulted in 33 chromosomes (Chr) (Figure S1), accounting for 99.98% of the total assembly (Table 4).In order to further determine the chromosome number of the Asian giant softshell turtle, about 1 mL blood was extracted from the neck blood vessels of four 1 to 3 years-old young individuals from the turtle domestication base of the Pearl River Fisheries Research Institute.After 72 h of culture, cells were collected for the preparation of chromosomes analysis as previously described 14 .A total of 43 metaphase mitotic cells with a higher number of chromosomes were obtained (Figure S2, Figure S3).Finally, we selected one of the chromosome division phases to make the standard karyotype map (Figure S4), and the karyotype formula was determined by the first eight division phase maps in Supplementary figure 1. Subsequently, Juicer and JucieBox v.1.8.8 were implemented to construct an interactive map and correct assembly errors visually, respectively.The heat map of the Hi-C assembly interaction bins displayed stronger signals around the diagonal than that of other positions.The phenomenon revealed that intensity of the interaction between adjacent sequences (diagonal position) was high, while the intensity of the interaction signal between non-adjacent sequences (non-diagonal position) was weak, which was consistent with the principle of Hi-C assisted genome assembly, proving a high quality and completeness of a genome assembly (Fig. 3).assessment of the genome assemblies.Two methods were used to check the completeness and quality of the assembly.First, BWA v0.7.10-r789 15 was used to compare the short sequences obtained by second-generation high-throughput Illumina sequencing with the reference genome, and the integrity of the assembled genome was evaluated by the alignment rate.Approximately 98.95% of clean reads were mapped to contigs, and 96.96% of clean reads were mapped to proper pairs (Table 5).Second, BUSCO v3.0 16 was used to search the 2,586 universal single-copy orthologous genes in vertebrata_OrthoDB v9 to evaluate the completeness, degree of fragmentation and missing genes of the genome assembly.Among the 2,586 prospective conserved core genes in the eukaryotic database, 2,478 (95.82%) and 54 (2.09%) were identified as complete and fragmented BUSCOs, respectively, indicating that the assembled genome had high integrity and validity, and could be used for further analysis (Fig. 4, Table 6).

Gene prediction and functional annotation.
With the help of LTR_FINDER 17 and RepeatScout 18 , a repeat sequence database of the genome was constructed, based on the principle of structure prediction and ab initio prediction.PASTEClassifier 19 was used to categorize different types of repetitive sequences.The ultimate repeat library was subsequently determined combined with Repbase 20 database.At last, RepeatMasker 21 software was conducted to predict the repeated sequences.We identified 993,573,618 bp sequences as repeats in the Repbase database, which covered 45.98% of the genome assembly (Table 7).The most common repetitive elements, RNA transposons (class I), account for 39.33% of the genome content, higher than DNA transposons (Class II) (Table 7).Long interspersed nuclear elements (LINEs) were the most abundant repeating element, followed by large retrotransposon derivatives (LARDs) and terminal inverted repeats (TIRs) (Table 7).In addition, several sequences classified as unknown repeats accounted for 2.48% of the genome assembly (Table 7).
A combination of three different strategies including ab initio prediction, homologous species prediction and transcriptome-based prediction were conducted to establish gene models.For gene predication based on ab initio, Genscan 22 , Augustus v2.4 23 , GlimmerHMM v3.0.4 24 , GeneID v1.4 25 , and SNAP version 2006-07-28 26 were employed to predict protein-coding genes.For homology-based prediction, GeMoMa v1.3.1 27,28 was used to align the assembled genome of P. cantorii with homologous species to predict the potential gene structures.Regarding RNA-Seq based prediction, Hisat v2.0.4 29     based on reference transcripts, and TransDecoder v2.0 and GeneMark-ST v5.1 31 were used for gene prediction.PASA v2.0.2 32 was used to predict Unigene sequences based on unreferenced assembly of transcriptome data.Finally, EVM v1.1.1 33 software was used to integrate the prediction results of the above three methods, and PASA v2.0.2 was used to modify the results (Table 8).
The predicted gene sequences were compared using BLAST V2.2.31 34 (-evalue1e-5) with functional databases such as NR 35 , KOG 36 , GO 37 , KEGG 38 and TrEMBL 39 .Gene functional annotation analysis, including KEGG pathway annotation analysis, KOG functional annotation analysis and GO functional annotation analysis, was performed.Finally, 21,833 protein-coding genes were successfully annotated.Of these predicted genes, 21,046 (~96.40%) were functionally annotated in at least one database, including GO, KEGG, KOG, TrEMBL and NR database (Table 9).Moreover, Blastn was used in the Rfam 40 database for whole-genome comparison to identify microRNAs and rRNAs.tRNAscan-SE 41 was used to identify tRNA.A total of 1595 tRNAs, 220 microR-NAs and 86rRNAs were identified.

technical Validation
Genomic integrity, fragmentation, and possible loss rates were measured using BUSCO V3.Among 2,586 prospective conserved core genes in the eukaryotic database, 2,478 (95.82%) and 54 (2.09%) were identified as complete BUSCOs and fragment BUSCOs, respectively, indicating that the assembled genome had high integrity and validity and could be used for further analysis (Fig. 4, Table 6).

Fig. 2 A
Fig. 2 A 21-mer distribution of the Illumina short reads to estimate genome size, ratio of repeat sequences and heterogeneity.The x-axis represents the sequencing depth of each unique 21-mer, and the y-axis represents the frequency of unique 21-mers.The k-mer depth value of the main peak is 56.After removing data with abnormal depth, a total of 131,112,165,086, 21-mers were used to further estimate the genome size; thus, we first estimated that the genome size of the Asian giant softshell turtle was 2.20 Gb.

Fig. 3
Fig. 3 Hi-C assembly of Chromosomal interaction maps.The 33 squares represent the constructed 33 chromosomes (Chr01 -Chr33) of the Asian giant softshell turtle.The color from light (low) to dark (high) indicates contact density of Hi-C interaction links.

Fig. 4
Fig. 4 Genome completeness of the Asian giant softshell turtle genomic sequence evaluated by BUSCO v3.0.The colours represent different kind of BUSCOs the Complete BUSCOs, Fragmented BUSCOs, Missing BUSCOs, Complete and duplicated BUSCOs and Complete and single-copy BUSCOs.

Table 1 .
The data statistic of k-mer analysis and heterozygosity in the Asian giant softshell turtle.

Table 2 .
Statistics of short-reads illumina in the Asian giant softshell turtle genomic survey.
8.1 was performed to filter

Table 3 .
Statistics of the genome assembly, Hi-C results and gene sent.

Table 4 .
Table summary of the assembled chromosomes of the Asian giant softshell turtle.

Table 5 .
The alignment of the Illumina reads to the Asian giant softshell turtle genome assembly.

Table 7 .
The detailed repetitive elements in the genome.

Table 8 .
Statistics of gene prediction results.

Table 9 .
Functional annotation from the genome assembly of the Asian giant softshell turtle.