De novo assembly and characterization of the Chinese three-keeled pond turtle (Mauremys reevesii) transcriptome: presence of longevity-related genes

Ethical approval

Procedures involving animals and their care were approved by the Animal Care and Use Committee of Anhui Normal University under approval number #20140111.

Sample collection and RNA extraction

Three adult female turtles were collected from our plant at Wuhu, Anhui, China. The liver, spleen, and skeletal muscle were collected and dissected. Tissue samples were stored at −80 °C. Total RNA was extracted from the tissues separately using a Trizol kit (Invitrogen, CA, USA) according to the manufacturer’s protocol. Extracted RNA was quantified with Nanodrop (Thermo, CA, USA) and the integrity and size distribution were checked with agarose gel electrophoresis. High-quality RNA from all three tissues was pooled for cDNA synthesis and sequencing.

cDNA library construction and sequencing

Two micrograms of pooled total RNA was used for cDNA library construction using TruSeq^® RNA LT Sample Prep Kit v2 (Illumina, CA, USA) according to the manufacturer’s protocol. We then prepared the synthesized cDNA for sequencing library construction by performing end-repair, 3′-end adenylation, as well as adapter-ligation and enrichment. Sequencing was performed using an Illumina Hiseq 2500 platform (Quail et al., 2008) at Genergy Bio-technology Co., Ltd. (Shanghai, China).

Sequence data processing and de novo assembly

The raw reads generated by the Illumina sequencer were saved as fastq format files. Adapter sequences were trimmed and low-quality reads were removed from the raw reads using Trim Galore (version 0.3.5) software. FastQC (version 0.10.1) was used to check the quality of pretreated data; reads that achieved a high Phred score (>28) were used for the assembly. We used Trinity (Grabherr et al., 2011), with a k value = 25, to perform de novo assembly on the trimmed and quality-checked reads. Sequence data were partitioned into many individual de Bruijn graphs, representing the transcriptional complexity at a given gene or locus. Each graph was processed independently to extract full-length splicing isoforms and to tease apart transcripts derived from paralogous genes. The final Trinity output was analyzed. Gene expression level was calculated using RSEM software (version 1.2.3) (Li & Dewey, 2011).

Blast against turtle’s reference proteomes

We used Blastx to query all M. reevesii transcripts to the proteomes of the green sea turtle, the western painted turtle, and the soft shell turtle, which were downloaded from the GenBank. The e-value cutoff was set to 1 × e⁻¹⁰ and the maximum target number was set to 20. The top results were selected as the annotation of the gene (termed as “unique protein”). A Venn diagram of the homologous genes across the three turtle proteomes was generated with VENNY (Oliveros, 2007).

Functional annotation

Sense and component strands of potential protein coding sequences (CDS) were predicted using Transdecoder in Trinity software v2.0.2 package (http://transdecoder.github.io/), based on a Markov model with default parameters. CDS were then translated into amino acid sequences with reference to a standard codon table. We used Blastp to search the Swissprot/Uniprot database (Balakrishnan et al., 2005) with our potential protein sequences as queries. The e-value cutoff was set to 1 × e⁻³. A gene name was assigned to each contig based on the top Blastp hit.

Gene ontology (GO) analyses (Ashburner et al., 2000) on all predicted protein sequences were conducted using InterProScan, set to default parameters (Zdobnov & Apweiler, 2001). The GO terms associated with each assembled sequence of the turtle transcriptome were then classified into biological processes, molecular functions, and cellular components. InterProScan was also used to predict the functional domains, signal peptides, and other protein characters by blasting to the Conserved Domain Database Interpro.

We then employed the KEGG Automatic Annotation Server (KAAS) with default settings (Moriya et al., 2007) to perform a KEGG pathway analysis (Kanehisa & Goto, 2000) on each contig. Telomere-related genes were filtered from the annotation results.

Analysis of candidate longevity-associated genes

Candidate genes associated with turtle longevity were screened from the logfile generated by Blastp. Two shelterin genes (i.e., pot1 and tin2) were not detected in this step. To search the lost but should be existed genes in M. reevesii transcriptome, we used released pot1 (XM_007063192.1) and tin2 (XM_005290151.2) genes to blastn against the primary assembled M. reevesii transcriptome data. And the pot1 and tin2 genes were finally screened out that escaped from functional annotation.

The filtered longevity-associated genes were then blasted to the NR database using Blastn, to verify the accuracy of previous annotations. Candidate gene fragments were aligned, overlapping fragments were assembled into one sequence, and separated fragments were concatenated into a single sequence. FPKM values of each fragment were also checked using the results from the RSEM software.

Screened genes were translated into amino acids, using MEGA 6.06 (Tamura et al., 2013) and according to a standard codon table. Amino acid sequences were then aligned with published data using the online alignment tool MAFFT (version 7) (http://mafft.cbrc.jp/alignment/software/) (Katoh et al., 2002), with default settings. BioEdit 7.2.3 (Hall, 1999) was used to display the alignment results.

Sequencing and de novo assembly

Total RNA extracted from liver, spleen, and skeletal muscle were used to generate a pooled cDNA sample and subsequently sequenced. A total of 160,998,396 paired-end raw reads of 100 bp length were generated and stored in fastq format (GenBank accession number: SRX1469958). We obtained 152,214,434 (94.5%) high-quality reads with an average length of 98.8 bp. The results of de novo assembly yielded 459,911 isoforms, which clustered into 230,085 transcripts with average length of 660 bp and median length of 342 bp. 190 transcripts encompassed over 100 isoforms and 6,051 transcripts encompassed over 10 isoforms with an average of two isoforms per transcript overall. Additionally, the highest number of isoforms found in one transcript was 1,048. These results may indicate the widespread existence of alternative splicing in M. reevesii. GC content for the entire final assembly and all protein coding sequences were 46.67% and 51.16%, respectively. Repetitive elements and microsatellites were also analyzed (File S1).

Comparison with turtle’s reference proteomes

The assembly quality of the M. reevesii transcriptome was assessed with a Blastx comparison to the reference proteomes of three turtles: the green sea turtle (Chelonia mydas), the Chinese soft-shell turtle (Pelodiscus sinensis), and the western painted turtle (Chrysemys picta bellii) (Fig. 1 and Table 1).x Chelonia mydas exhibited the highest degree of similarity. The comparison to P. sinensis and C. picta yielded slightly more unique proteins than transcripts hits, but overall, we obtained a total of 22,776 positive transcripts hits.

Transcripts expression level

The expression of each transcript was quantified using RSEM software, set to default parameters. Statistical results are shown in Table 2. The transcripts with the highest levels of expression were related to metabolism and translation activity.

Functional annotation

The results of our functional annotation revealed 42,918 (18.7%) transcripts that were predicted to potentially code for proteins. These were then translated into amino acid sequences by referring to a standard codon table. After a blastp search on the resultant amino acid sequences, we obtained a total of 534,077 hits, the best of which corresponded to 18,846 unique protein accessions in the Swissprot/Uniprot database. Length distribution for all opening reading frames ranged from 40 bp to 34,967 bp, with an average length of 874 bp (Fig. 2). The Blastp top-hit species distribution of gene annotations showed the highest homology to Chelonia mydas (6,654 annotation results) and Pelodiscus sinensis (5,747 annotation results) (Fig. S1). Those two species supplied 65.8% of all annotation information due to their close phylogenetic relationships to M. reevesii.

The results of a gene ontology (GO) analysis assigned 11,695 unique proteins to 4,031 terms for biological processes, molecular functions, and cellular components. Within biological processes, metabolic (26%) and cellular processes (26%) were the most well-represented. Next, the majority of the proteins assigned to molecular functions were associated with binding (85%). Finally, within cellular components, cell (33%) and membrane (31%) proteins were the most well-represented (Fig. 3).

The most abundant conserved protein domain found in our data was the zinc finger C2H2 domain, followed by the immunoglobin-like fold and fibronectin-III. Zinc finger-associated conserved protein domains are 8% of all conserved domains, and 62.7% of those are associated with zinc finger C2H2 (Table 3).

The KEGG pathway analysis annotated 4,486 transcripts into 338 pathways, with 3.5 KEGG pathways per transcript on average. Of all the annotated sequences, a large portion (1,581, or 35.2% of 4,486) were related to metabolism, specifically of carbohydrates (332 sequences), lipids (247 sequences), and amino acids (291 sequences). Next, 855 sequences (19.1%) were involved in signal transduction. The most well-represented was the PI3K-Akt signaling pathway (ko04151; 224 sequences), which could be activated by IGF-1, followed by the MAPK signaling pathway (ko04010; 159 sequences). Finally, 483 sequences (10.8%) were associated with the immune system, including the T cell (ko04660) and B cell (ko04662) receptor signaling pathways.

Longevity-related genes

Candidate genes related to longevity were screened out from the total pool of annotated genes. Specifically, fragments of tert, tep1 and six (trf1, trf2, tpp1, pot1, tin2 and rap1) shelterin proteins were found (File S2). Genes were translated referring to standard codon table and aligned with public database.

Tert gene, encoding the catalytic subunit of telomerase enzyme, has been drawn a lot of attention by biologists. We obtained a 3,817 bp of tert gene fragment including two conserve domains in 3^′′ terminal: RNA binding domain of telomerase (TRBD) and reverse transcriptase (RT) domain are essential for the gene function. The tert gene sequence length and identity were compared with other turtles and human (Table 4). Nucleotide conservation is interestingly higher than amino acid conservation, the similar results were also found in Nothobranchius furzeri (Hartmann et al., 2009).

FPKM values are approximately 10 of all aforementioned genes and genes that were closely related in function (File S3). The similar expression levels indicated that the candidate genes were accurately identified.

Profile of the M. reevesii transcriptome

More than 15-Gb high-quality data were generated with Illumina sequencing in this study. A total length of 16,478,555 bp (0.75% of 2.2 Gb) for 42,918 assembled transcripts were annotated that encoded for mRNA in M. reevesii transcriptome, functional annotation of the mRNAs provided fully information of transcripts and whole profile of the transcriptome that could be used in further study. The most abundant domain, C2H2 zinc finger, was consistent with the results found in mammalians and humans, the fact that the motif is the most prevalent and the largest sequence-specific DNA-binding protein family (Lander et al., 2001; Tupler, Perini & Green, 2001). The conservation and evolution of this domain could be further studied (Englbrecht, Schoof & Böhm, 2004).

Lastly, our functional annotation results also revealed a small amount of contaminated transcripts associated with viruses. While such data are negligible, they should nonetheless be excluded in future studies.

The active expression of telomere-associated genes may contribute to turtle longevity

Telomere shortening is now considered the molecular clock that triggers cell and organismal senescence (Harley & Goldstein, 1978). To prevent premature shortening, telomere length is maintained by telomerase, previously considered to be inactive in human somatic cells (Kim et al., 1994; Shay & Bacchetti, 1997). Using RNA-seq technology, we were able to find two coding genes of the telomerase complex—TERT and TEP1—in the M. reevesii transcriptome. Both genes are essential for proper telomerase function: transient expression of TERT has been found to reconstitutes telomerase activity (Weinrich et al., 1997), while TEP1 is a component of the ribonucleoprotein complex (Poderycki et al., 2005). While TEP1 has been closely linked to telomerase activity (Nakayama et al., 1997), other research has thrown doubt on this connection (Uchida et al., 1999). Thus, the exact function of TEP1 is a subject for further research. The presence of these two genes might implies that activated telomerase is present in turtle somatic tissues to protect the turtle from senescence.

We were able to detect all the reported six genes that encode shelterin proteins in our data. Expression of these genes could prevent the telomeres from shortening during cell division and thus were predicted having a significant effect to the longevity of this species.