SMRT sequencing of the full-length transcriptome of the Rhynchophorus ferrugineus (Coleoptera: Curculionidae)

Samples selection and preparation

All R. ferrugineus samples used in this study were collected from the Coconut Research Institute, Chinese Academy of Tropical Agricultural Sciences, Wenchang Hainan, China. Samples were divided into larva, pupa, female adult and male adult. The whole body was collected from 12 R. ferrugineuss (3 from 7th instar larvae, 3 from Pupae, 3 from female adults, 3 form male adults). All samples were harvested and frozen in liquid nitrogen and stored at −80°C for further experiments.

RNA extraction and SMRT sequencing

Total RNA samples were isolated using the RNeasy Plus Mini Kit (Qiagen, Valencia, CA, USA). Then 1% agarose gels was used to detect RNA degradation and contamination, and Nanodrop (NanoDrop products, USA) was used to check RNA purity (OD 260/280). RNA concentration and integrity were accurately evaluated using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA) and Agilent 2100 (Agilent Technologies, USA), respectively. For PacBio Iso-Seq, the total RNA samples from three developmental stages (larva, Pupa, female adult and male adult) were mixed together for the following experiments. A total of 3 µg of mixed RNA was sequenced on the Pacbio Sequel platform (Pacific Biosciences, CA, USA) in accordance with the manufacturer’s instructions. Then, according to the Isoform sequencing protocol (Iso-Seq), the Iso-Seq library was prepared by the Clontech SMARTer PCR cDNA synthesis kit (Clontech, CA, USA) and the BluePippin size selection system protocol described by Pacific Biosciences (PN 100-092-800-03). For Illumina RNA-Seq, twelve libraries of three developmental stages (larva, pupa, male adult and female adult) RNA samples were prepared and sequenced respectively. For each sample, a total amount of 1.5 µg RNA was used for short reads sequencing on Novaseq 6000 platform. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendations and index codes were added to attribute sequences to each sample. The sequencing work reported in this work was performed by Novogene technology co. (Beijing, China). PacBio Iso-Seq and Illumina RNA-seq data generated from R. ferrugineus are available from the NCBI SRA database under project number PRJNA598560.

Data processing and error correction of PacBio Iso-Seq reads

Firstly, SMRTlink 6.0 software was used to process the sequence data. Immediately after, the cyclic consensus sequence (CCS) was generated from the subread BAM files (parameters: min_length 50, max_drop_fraction 0.8, no_polish TRUE, min_zscore-9999.0, min_passes2, min_predicted_accuracy 0.8, max_length 15000) and a CCS.BAM file was output. The generated BAM files were divided into full-length and non-full-length reads using pbclassify. Finally, input the full and non-full length fasta files into the clustering step, which performs the isoform level clustering, and then uses Quiver (parameters: hq_quiver_min_accuracy 0.99, bin_by_primer false, bin_size_kb 1, qv_trim_5p 100, qv_trim_3p 30) for the final arrow polishing. Full-length transcripts were corrected using Illumina RNA-seq data with the software LoRDEC (Salmela & Rivals, 2014). The redundancies in the corrected transcript were then removed using the CD-HIT-EST (parameters: -c 0.95 -T 6 -G 0 -aL 0.00 -aS 0.99) program to obtain the final transcript for subsequent analysis (Fu et al., 2012).

Functional annotation of transcripts

Transcripts function was annotated based on the following databases:NR (NCBI non-redundant protein sequences) (Deng et al., 2006), NT (NCBI non-redundant nucleotide sequences), Pfam (Protein family) (Finn et al., 2014), KOG (Clusters of Orthologous Groups of proteins) (Koonin et al., 2004), Swiss-Prot (A manually annotated and reviewed protein sequence database) (Apweiler et al., 2004), KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa et al., 2004) and GO (Gene Ontology) (Ashburner et al., 2000). The BLAST software with E-value ≤1e−5 was used for NT database analysis. The Diamond v0.8.36 software with E-value ≤1e−5 was analyzed in NR, KOG, Swiss-Prot and KEGG annotations. The Hmmscan procedure was used in the Pfam database, and GO function categories were performed by Blast2GO v2.5 based on Pfam annotation.

CDS prediction and SSR analysis

The ANGEL pipeline is a long-read implementation of ANGLE that is used to determine protein coding sequences from cDNAs. We use the confidence protein sequences of R. ferrugineus or closely related species for ANGLE training, and then run the ANGLE predictions for given sequences (Shimizu, Adachi & Muraoka, 2006). At the same time, the MISR (http://pgrc.ipk-gatersleben.de/misa/misa.html) was used to identify and locate the SSR of the transcriptome.

Identification of TFs and lncRNAs

Animal transcription factors were performed by the animal TFDB 2.0 database (Zhang et al., 2015). For species included in the database, if they were not Ensembl geneid genes, they would be screened by BLASTX with the known transcription factors protein sequence of the species in the database, and if they were Ensembl geneid, they would be screened directly. For species not included in the database, hmmsearch was used to identify them according to pfam files of the transcription factor family. LncRNA of the transcriptome were predicted by using Coding-Non-Coding-Index (CNCI) (Altschul et al., 1997), Coding Potential Calculator (CPC) (Kong et al., 2007), Pfam-scan (Finn et al., 2016) and PLEK (Li, Zhang & Zhou, 2014) to predict the coding potential of transcripts. Firstly, PLEK SVM classifier with default parameters of—minlength 200 and CNCI with default settings were performed to evaluate coding potential; Secondly, CPC and NCBI eukaryotic protein database were used for BLAST comparison (E-value < 1e−10 setting) to search transcripts; Finally, homologous search of hmmscan was performed for the transcription sequences predicted by the three software with the Pfam database, and the protein family domains were recorded, with the default parameter of -E 0.001-domE 0.001.

AS analysis

To obtain alternative splicing (AS) events for R. ferrugineus, transcripts were further processed using Coding GENome reconstruction tool (Cogent v3.1, https://github.com/Magdoll/Cogent). Generally, Cogent first divides the input fasta file into chunk_size blocks, and then calculates the k-mer configuration file. Then, the De Bruijn graph was used to further reconstruct each transcription family into one or more unique transcription models (called UniTransModels). Finally, gmap-2017-06-20 was conducted to map the adjusted transcripts to UniTransModels. Splicing junctions detection was performed on transcripts mapped to the same UniTransModels, and transcripts with the same splice junctions were collapsed. Meanwhile, the transcriptional isoforms of UniTransModels have collapsed transcripts with different splicing junctions. AS events were detected with SUPPA (https://github.com/comprna/SUPPA) using default settings.

The full-length sequences of R. ferrugineus using PacBio sequencing

The full-length transcriptome of R. ferrugineus was generated using the PacBio Sequel platform on the pooled RNA of twelve R. ferrugineus samples. The results showed that PacBio Sequel platform generated a total of 454,369 circular consensus sequences (CCSs) with a full length reads of 362,466. The nonfull-length (nFL) sequences was 81,424 and the full-length non chimera (FLNC) reads number was 330,973 with an average length of 2,332 bp. Twelve samples were sequenced by Illumina Novaseq 6000 respectively, and a total of 642,179,304 raw reads and 625,983,256 clean reads (97.48%, 93.91 G) were obtained. PacBio Sequel platform produced a total of 10,172,136 subreads and 181,405 consensus reads (16.67 G, with an average length of 2,305 bp, an N90 of 1,327 bp and an N50 of 2,790 bp), which were then corrected using the Illumina reads (after correction, the average length was 2,302 bp, N90 was 1,321 bp, and N50 was 2,785 bp) and ubsequently removing redundancy via the CD-Hit program, the consensus transcripts were finally clustered into a total of 63,801 transcripts for subsequent analysis. We found that the main length distribution range of unigenes was 0.5–6 k (Fig. 1).

Gene annotation of R. ferrugineus

To obtain a comprehensive functional annotation of R. ferrugineus transcriptome, we annotated 63,801 transcripts with seven databases, including Swiss-Prot, KOG, GO, NR, NT, Pfam, and KEGG. In total, 50,280, 40,109, 47,197, 33,511, 27,707, 27,253 and 27,707 transcripts were annotated in the NR, Swiss-Prot, KEGG, KOG, GO, NT and Pfam databases, respectively. Moreover, 54,999 transcripts were annotated in at least one database and 12,508 transcripts were annotated in all databases (Fig. 2).

NR is a non-redundant protein database characterized by its comprehensive content and the inclusion of species information in the annotated results, which can be used for the classification of homologous species. Aligned all transcripts in the NR database, the results showed that 50,280 transcripts were annotated in the NR database and the top five most annotated in NR database were Dendroctonus ponderosae, Anoplophora glabripennis, Bactrocera tryoni, Tribolium castaneum and Aethina tumida (Fig. 3).

KOG is a protein database created and maintained by NCBI, which is constructed according to the phylogenetic relationship of coding proteins in complete genomes of bacteria, algae and eukaryotes, and is widely used to predict the function of sequences. The KOG functional classification of R. ferrugineus transcripts was shown in Fig. 4. The results indicated that a total of 33,511 genes categorizing into 26 categories were annotated in KOG database. The first six largest groups among these categories were transcription, general function prediction only, function unknown, signal transduction mechanisms, cytoskeleton and posttranslational modification (protein turnover and chaperones), respectively.

Full-length transcripts of red palm weevil were annotated with GO database, and 27,707 annotated transcripts were successfully divided into three categories: biological processes, molecular functions, and cellular Components (Fig. 5). In biological process, the cell process accounts for the largest proportion, followed by metabolic process and single-organism process. In addition, we also found that part of the genes was annotated into biological regulation, regulation of biological process, localization, response to stimulus and signaling terms. In cellular component, the genes involved in cell, cell part, organelle, membrane, membrane part and macromolecular complex were the most. In molecular function categories, binding, catalytic activity and transporter activity were identified as the most abundant terms.

KEGG Pathway analysis can be used to systematically analyze the metabolic pathways of gene products and compounds in cells and the functions of these gene products. In the KEGG classification of R. ferrugineus, human diseases, metabolism and organismal systems were the top three categories with higher proportions (Fig. 6). Briefly, a total of 11,950 genes were involved in human disease related pathways, in which 1,668 genes were predicted to infectious disease: viral, 1,396 genes were predicted tonfectious diseases: bacterial and 2,777 genes were predicted to cancers: overview. A total of 10,294 of the annotated genes were classified as belonging to organismal systems related pathways, in which the nervous system (1,016 genes), immune system (1,984 genes), digestive system (993 genes) and endocrine system (2,500 genes) were the top four pathways with the most abundant genes. In addition, a total of 7,979 annotated genes were involved in the Metabolism pathway. The most abundant pathways were carbohydrate metabolism (1,478 genes) and amino acid metabolism (1,180 genes). Regarding Environmental Information processing, genes were involved in the signal transduction (4,741 genes), signaling molecules and interaction (304 genes) and membrane transport (196 genes). At the same time, a lower number of genes are annotated to Cellular Processes and Genetic Information Processing.

CDS prediction

The CDS is a sequence encoding a protein product that corresponds exactly to the codon of a protein. In the sequencing results of the full-length transcriptome, predicting the protein coding region contributes to the preliminary analysis of the gene and is also the basis for subsequent protein structure analysis. For red palm weevil, ANGEL software was performed to carry out CDS prediction analysis on the obtained full-length transcriptome sequence, and the results showed that the main distribution range of CDS length was 0 ∼2,500 nt (Fig. 7).

Transcription factors identification

Transcription factors are an important part of the transcriptional regulatory system. Using the present data of R. ferrugineus, 2,084 transcription factors were predicted, and Zf-C2H2 (570,27.35%), ZBTB (476,22.84%), TF_bzip (101,4.85%) and bHLH (85, 4.08%) were the top four transcription factor families (Fig. 8). These transcription factors will lay the foundation for exploring the role of the regulatory mechanism of red palm weevil.

SSR discovery

Simple sequence repeats is a group of repeated sequences consisting of several nucleotides (1∼6) with repeat units up to dozens of nucleotides. The repeats are short in length and widely distributed uniformly in eukaryotic genomes. In our analysis, MISA software (version 1.0, default parameter) was applied for SSR detection of transcriptome. The results showed that the minimum number of repetitions of each unit size is 1–10, 2–6, 3–5, 4–5, 5–5, 6–5. In total 66,230 SSR loci were identified in this transcriptome; mono nucleotide motifs (49,898, 75.34%) were the most abundant type of SSR locus, followed by di nucleotide motifs (12,662, 19.12%), tri nucleotides(3,377, 5.09%), tetra nucleotides (192,0.29%), penta nucleotides (33, 0.05%) and hexa-nucleotide motifs (68, 0.10%) (Fig. 9).

LncRNA prediction

LncRNA (long-chain noncoding RNA) is a class of RNA molecules whose transcripts are more than 200 nt in length and do not encode proteins. For R. ferrugineus, the numbers of lncRNAs identified from transcriptome by CNCI, CPC, PLEK and Pfam were 32,552, 17,481, 18,897 and 34,066, respectively (Fig. 10). The intersection of these four results produced 9,618 lncRNA transcripts. Meanwhile, the length distribution density of mRNA was compared with the predicted lncRNA (Fig. 11).

Alternative splicing analysis from full-length transcriptomes

Since R. ferrugineus had no reference genome, we used Cogent (Coding genome reconstruction tool) to reconstruct genes using high-quality full-length transcriptome sequences to generate UniTransModels. UniTransModels were used as a reference sequence to describe the types of AS events and the number of corresponding genes. The results indicated that a total of 2,184 UniTransModels-based AS events in R. ferrugineus were detected. Briefly, six main AS events (alternative 3′ splice sites, Mutually exclusive exons, Skipping Exon, alternative 5′ splice sites, Retained introns and Alternative First Exons) were identified. Retained introns (RI) were identified as the most abundant event, accounting for 6.14% (134) of all events. The other five types of AS events account for less than 2% of all detectable events. The number of two kinds of events, alternative 3′ splice sites (32,1.47%) and alternative 5′ splice sites (34,1.56%), were slightly higher than those of skipped exons (7,0.32%), mutually exclusive exons (1,0.05%) and Alternative first exons (20,0.92%).