Transcriptome analysis between invasive Pomacea canaliculata and indigenous Cipangopaludina cahayensis reveals genomic divergence and diagnostic microsatellite/SSR markers

Sequencing data and de novo assembly

Using Illumina paired-end sequencing technology, the transcriptome sequencing produced 65,198,546 reads with a total length of 6.5 Gb for C. cahayensis, which generated 161,941 contigs and 151,518 unigenes (Table 1). For P. canaliculata, 94,808,488 reads were obtained, and 94,518 contigs and 76,082 unigenes were generated (Table 1). Using the SOAP de novo assembly program, high quality reads were assembled into 160,256 contigs longer than 200 bp, with a mean length of 1,080 bp and a N50 of 1,004 bp for the native C. cahayensis. For P. canaliculata, 94,518 contigs longer than 200 bp, with a mean length of 916 bp and a N50 of 1,854 bp were generated. In C. cahayensis, the lengths of 104,713 (65.34%) of the contigs ranged from 200 to 500 bp, 28,918 (18.04%) contigs ranged from 500 to 1,000 bp, and 15,191 (9.50%) contigs ranged from 1000 to 2,000 bp; the remaining contigs were longer than 2,000 bp (Figure 1). In P. canaliculata, the lengths of 41,544 (43.95%) of the contigs ranged from 200 to 500 bp, 19,289 (20.41%) contigs ranged from 500 to 1000 bp, and 17,619 (18.16%) contigs ranged from 1000 to 2,000 bp; the remaining contigs were longer than 2,000 bp. The related data were submitted to the NCBI data under accession numbers: SRA191276 (P. canaliculata) and SRA192725 (C. cahayensis).

Figure 1

Assessment of transcriptome assembly quality of Cipangopaludina cahayensis (A) and Pomacea canaliculata (B).

Functional annotation

To annotate the C. cahayensis and P. canaliculata sequences, searches were conducted against the NCBI non-redundant protein (Nr) database, the Swiss-Prot protein database, Cluster of Orthologous Groups (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) database using BLASTX (E-value ≤ 1 × 10⁻⁵). The alignment results were used to predict unigene transcriptional orientations and coding regions. Gene ontology (GO) terms were assigned to the annotated sequences and 14,864 sequences from C. cahayensis and 56,300 sequences from P. canaliculata were categorized into the three GO categories, biological process, cellular component, and molecular function (Figure 2). We found that the distribution and percentages of the assigned gene functions were similar in both species. In the biological process category, death (22.1%) was prominent, while in the molecular function category, cell (30%–31%) and cell part (30%–31%) were prominently represented. In the cellular component category, binding (47.8%–49%) was predominant, followed by catalytic activity (36%). Overall, the transcriptome sequencing yielded a great number of unique genes in the two species, in agreement with similar results reported in other species [20]. Several differences were noted between the two species, with more genes noted in P. canaliculata (56,300 genes) compared with in C. cahayensis (14,864 genes). Furthermore, the percentage of genes annotated as metabolic process/pigmentation under the biological process category was higher in P. canaliculata (15.7%/7.46%) compared with C. cahayensis (7.93%/1.6%), implying a possible relation to various environmental stressors. Moreover, the percentage of genes annotated as metallochaperone activity and translation regulator activity under the cellular component category was much higher in P. canaliculata compared with C. cahayensis. These results indicated that P. canaliculata might contain additional genes that are able to confer high competitiveness or strong resistance to envrionmental stress compared to C. cahayensis.

Figure 2

Comparing functional annotations of contigs between Cipangopaludina cahayensis (red) and invasive Pomacea canaliculata (blue) transcriptome. The distribution of gene ontology (GO) terms is given for each of each of the three main GO categories (biological process, molecular function, and cellular component).

Furthermore, all of the C. cahayensis and P. canaliculata unigenes were subjected to functional prediction and classification using the COG database. The unigenes were assigned to 25 COG categories (Figure 3), among which “general function prediction” represented the largest group (4,081 (17.9%) genes for C. cahayensis; 4,346 (19%) genes for P. canaliculata). For C. cahayensis, the next most represented category was translation, ribosomal structure and biogenesis (1915 (8.41%) genes), while for P. canaliculata, replication, recombination and repair (1,883 (8.23%) genes,) was the next most represented category.

Figure 3

Clusters of orthologous group (COG) classifications for Cipangopaludina cahayensis (A) and Pomacea canaliculata (B) transcriptome. All unigenes were aligned to COG database to predict and classify possible functions.

To identify differentially regulated biological pathways between C. cahayensis and P. canaliculata, the annotated unigenes were mapped to reference pathways in the KEGG database [31]. We found that 13,351 C. cahayensis unigenes mapped to 276 pathways and 13,808 P. canaliculata genes mapped to 240 pathways, with different pathway associations between the two species. In C. cahayensis, the largest number of genes included cancer (577 (4.32%) genes; pathway: ko05200), focal adhesion (496 (3.72%) genes; pathway: ko04510), ubiquitin mediated proteolysis (427 (3.2%) genes; pathway: ko04120), and Huntington’s disease (333 (2.49%) genes; pathway: ko05016). In P. canaliculata, the predominant pathways were metabolic (2241 (16.23%) genes; pathway: ko01100), cancer (530 (3.84%) genes; pathway: ko05200), focal adhesion (415 (3.01%) genes; pathway: ko04510) and Huntington’s disease (348 (2.52%) genes; pathway: ko05016). Collectively, these transcriptome sequences and pathway annotations provide an essential resource for further screening and expression analysis of candidate genes related to the invasive abilities of P. canaliculata.

Analysis of protein families and genes

A total of 15,632 protein families were identified based on sequence similarities (Figure 4); 13,490 families for C. cahayensis and 13,453 families for P. canaliculata. When the transcriptomes of the two species were compared, a total of 11,312 protein families were found to be conserved between the C. cahayensis and P. canaliculata transcriptomes, and 2142 and 2178 families for P. canaliculata and C. cahayensis, respectively, were found to be differentially expressed. Some of the differentially expressed proteins may be responsible for the unique features of each of these species. An enriched analysis of the GO terms assigned to the 11,312 conserved protein families, identified 12 protein families that were significantly enriched (Table 2), including RNA transport (380 (2.6%) genes), spliceosome (383 (2.62%) genes), and endoplasmic reticulum protein processing (358 (2.45%) genes), which are related to protein transportation and metabolism. The finding that GO terms related to protein transportation and metabolism were enriched is inconsistent with the results reported for other invasive species such as Bemisia tabaci [32], possibly suggesting the critical roles of these pathways in these two species. We identified a total of 12 protein families (p-value < 10⁻⁵) encoded by the differentially expressed genes in C. cahayensis (Table 3), including those assigned to pathways pertaining to cancer (97 (6.92%) genes), toxoplasmosis (87 (6.21%) genes), and apoptosis (71 (5.06%) genes). In P. canaliculata, we identified a total of 17 protein families (p-value < 10⁻⁵) encoded by the differentially expressed genes, including those assigned to pathways pertaining to lysosomes (84 (4.02%) genes), vitamin digestion and absorption (71 (3.4%) genes), ECM-receptor interaction (57 (2.73%) genes), and metabolism of xenobiotics by cytochrome P450 (49 (2.35%) genes). We used reads per kilobase per million mapped reads (RPKM) to analyze the expression levels of P. canaliculata genes and identified 20 annotated genes with very high expression levels (RPKM > 2000), which were predicted to be involved in cell and protein structure (ferritin [Swiss-Prot: C7TNT3] and augerpeptide hhe53 [Swiss-Prot: P0CI21]) and ribosomes (60S ribosomal proteins and 40S ribosomal protein S8) (Table 4).

Figure 4

Protein families from the transcriptomes of Cipangopaludina cahayensis and Pomacea canaliculata. Protein families were identified for all the translated genes of the two transcriptomes using Blastp and a Markov Cluster algorithm (MCL), with the total number of protein families belonging to each category listed in the figure for the 11,312 protein families belonging to the two transcriptomes.

P. canaliculata has become an important pest in China and has exhibitied resistance to pesticides such as metaldehyde and niclosamide ethanolamine salt [33-35]; however, the molecular mechanisms underlying this resistance are still unclear. To detect unique resistance-related sequences, the unigenes were edited manually to remove redundant and overlying short sequences and the edited sequences were then used to identify genes encoding proteins related to the metabolism of pesticides. We identified P450 cytochromes (CYPs), a major family of enzymes involved in detoxification and metabolism, as potential major detoxification component proteins [36-38]. Previous studies have reported a correlation between increased exposure to metabolic neurotoxic pesticides and over-expression of P450 genes in many pest species [39-46]. In our study, 210 P450-related sequences were identified in P. canaliculata and only 159 were found in C. cahayensis, indicating that the number of P450 genes may be one of the contributory factors to pesticides resistance in P. canaliculata. While the number of P450 genes detected is not necessarily related to gene expression levels, an increased gene number of genes may increase metabolic enzyme detoxification activity, and contribute to the development of a progressive resistance in P. canaliculata. These findings will enhance the understanding of pesticide metabolism and help in the development of effective treatments for invasive species. To investigate the relationship between the P450 sequences from both species a phylogenetic tree was constructed using the neighbor joining (NJ) method in conjunction with bit-score values. Sixty of the sequences showed high homology and were classified into the CYP2, CYP3, and CYP4 families based on their similarity to sequences in the Nr database. These sequences clustered into three clades in the phylogentic tree that corresponded to the same three P450 families (Figure 5). We found a high concentration of P. canaliculata genes in the CYP4 family, possibly implying that these genes played important roles in the metabolism of pesticides in this invasive species. While these finding are insightful, they need to be examined further using RACE technology and RT-PCR before they can be accepted.

Figure 5

Neighbor-joining phylogenetic analysis of cytochrome P450 from Cipangopaludina cahayensis (CC) and Pomacea canaliculata (PC). CYP represent cytochrome P450.

Detection of intraspecific genetic variation

EST-SSRs serve as effective molecular markers for genetic mapping, comparative genomics and population genetic analysis in many invasive species. Characterization of EST-SSRs may enable breakthroughs in the detection of cryptic species, aid in defining the number and location of establishment events, and help trace the routes of alien species as they spread into new regions [47-51]. Compared with traditional methods, EST-SSRs are more transferable and advantageous than random genomic SSRs, enabling improved genetic studies related to population genetics [52]. Unitl now, only a few SSRs have been identified in P. canaliculata [20,53], which has hampered marker applications in this species. To further understand the invasive and adaptive mechanism in P. canaliculata, six P. canaliculata samples were collected from three. invasive regions/habitats in mainland China and examined for polymorphisms. A total of 16,717 potential SSRs were identified. As shown in Table 5, the di-nucleotide repeats were the most abundant (10,554, 63.1%), followed by tri-(4,480, 26.8%), tetra- (1,021, 6.10%), hexa-(341, 2.0%), and penta-nucleotide (321, 1.9%) repeats. The most abundant repeat combination was AG/CT (40.4%), followed by AT/AT (18.3%), AAG/CTT (7.8%), AAT/ATT (4.7%), AC/GT (4.0%) and ATC/ATG (3.4%) (Figure 6A). Based on the SSR-containing sequences, 8,428 SSR primers were developed and 100 SSRs (Additional file 1: Table S1) were selected to design EST-SSR primers based on the information (name and longer length of gene identified). Of the 100 SSRs examined by PCR amplification, 26 (26.0%) PCR products exhibited more than one band, which may have resulted from high heterozygosity, while the others SSRs generated bands of the expected length. In total, 143 amplicons were detected from the 100 primer pairs. The number of amplicons per primer pair ranged from one to three, with an average of 1.43 (Figure 6B). To estimate EST-SSR marker novelty, the amplicons were evaluated against previously reported P. canaliculata markers [20,53]. We found that the 100 EST-SSR markers had not been reported previously. Thus, other EST-SSR primers can be designed from the 8,428 identified EST-SSR to contribute further to the characterization of the invasive and adaptive processes. P. canaliculata and C. cahayensis have very similar morphological features, especially at the immature stages, which makes early identification difficult. Therefore, a molecular means for the identification and characterization of these two species is essential. Using the P. canaliculata SSR primers, we identified a unique amplicon (FSLssr64; Additional file 1: Table S1) that was present in P. canaliculata but absent in C. cahayensis (Figure 6C). Thus, FSLssr64 could serve as a species-specific molecular marker to distinguish these two species and aid in the prevention and detection of invasive P. canaliculata in different regions.

Figure 6

Frequencies and polymorphisms of classified SSR repeat types and molecular characterization of Pomacea canaliculata. (A): The graph shows the frequency of each repeat motif classified, considering the sum of the frequencies for complementary sequences (for example, the sum of frequencies for the dinucleotides AC and its complementary GT). (B) Polymorphism and validation of a subset of the microsatellite primer pairs for six P. canaliculata samples by agarose-gel profiling. 1–6 represent GZ1, GZ2, HN1, HN2, SG1, and SG2, respectively. (C) The SSR primer (FSLssr64) for species-specific identification between P. canaliculata and C. cahayensis.

Ethics statement

This study was approved by the Animal Care and Use committee of Aquatic Invasive Risk Assessment Center, Pearl River Fisheries Research Institute, Chinese Academy of Fishery Sciences.

Sample collection, RNA extraction, and next generation sequencing

P. canaliculata (20–25 mm shell length; 25.23 ± 0.34 g; 10 individuals) and C.cahayensis (20.4–23.2 mm shell length; 22.43 ± 0.46 g; 10 individuals) were collected without the use of chemicals and grown in the Aquatic Invasive Risk Assessment Center, Pearl River Fisheries Research Institute, Chinese Academy of Fishery Sciences, Guangzhou, China. Tissues samples from the foot, muscle, liver, and kidney were rinsed separately with water pretreated by diethyl pyrocarbonate to cleanse the samples and inactivate RNases [32]. Total RNA of each sample was extracted using a Trizol Kit (Promega) according to the manufacturer’s instructions. RNA quality was assessed using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and RNase-free agarose gel electrophoresis, with the total RNA concentration measured using a 2100 Bioanalyzer. Equal amounts of RNA from each sampled tissue were combined for subsequent experiments and RNA purity was assessed at absorbance ratios of OD_260/280 and OD_260/230. RNA integrity was confirmed by 1% agarose gel electrophoresis.

De novo assembly and gene annotation of Illumina reads

Transcriptome de novo assembly was carried out with the short-read assembly program Trinity [54]. The Trinity program has three independent modules: Inchworm, Chrysalis, and Butterfly. Inchworm assembled the RNA sequencing data into unique transcripts that we called Inchworm contigs; Chrysalis clustered the Inchworm contigs, then constructed complete de Bruijn graphs for each cluster and partitioned the full read set among these disjoint graphs; and Butterfly processed the individual graphs in parallel, tracing the paths based on reads and pair-end information, ultimately reporting full-length transcripts for alternatively spliced isoforms. After assembly, the TIGR Gene Indices clustering tools (TGICL) [55] were used to cluster and remove redundant transcripts. The remaining sequences after TGICL clustering were defined as unigenes. BLASTX searches (E-value < 10⁻⁵) were conducted to screen the unigenes against the Nr database (http://www.ncbi.nlm.nih.gov/), Swiss-Prot protein database (http://www.expasy. ch/sprot/), the KEGG pathway database (http://www.genome.jp/kegg/), and COG database (http://www.ncbi.nlm.nih.gov/cog/). High scoring alignments were used to determine the unigene sequence direction. When alignment results varied between databases, the results from the Nr database were preferentially selected, followed by the Swiss-Prot, KEGG and COG databases. When a unigene sequence did not match entries in any of these databases, ESTScan was used to predict the coding regions and determine sequence directionality.

Functional annotation and differential expression analysis of unigenes

Unigene sequences were aligned to the protein databases (listed above) using BLASTX (E-value < 10⁻⁵) and to the nucleotide sequence database Nt (E-value < 10⁻⁵) using BLASTN to obtain both protein and functional annotation information. Based on the annotations in the protein databases, Blast2GO [56] was used to obtain GO annotations for the aligned unigene sequences and the Web Gene Ontology Annotation Plot (WEGO) software [57] was used to establish GO functional classifications for all unigenes. The unigenes were aligned to the COG database to predict and classify possible functions and the KEGG database was used [31] to obtain pathway annotations (E-value threshold 10⁻⁵). RPKM was used to calculate unigene expression levels, which eliminated the influence of gene length and sequencing level on the estimation of gene expression

Analysis of genes related to pesticide

BLASTX searches against the Nr database (E-value < 10⁻⁵) were used to detect genes related to pesticide resistance. Sequences that returned redundant BLAST results or showed high sequence homology were eliminated and presumed to be allelic variants or different parts of the same gene. Thirty P450 gene sequences (Additional file 1: Table S1) with a range of bit-score values were identified and aligned using MUSCLE [58], and their phylogenetic relationships and genotype classifications were determined using MEGA 5 software [59]. The NJ method [60] was used to create phylogenetic trees, with positions containing alignment gaps or missing data eliminated via pairwise deletion. Tree branch strength was evaluated via a bootstrap analysis of 1000 replication trees.

Development and detection of EST-SSR markers

MIcroSAtellite (MISA) (http://pgrc.ipk-gate-rsleben.de/misa/) was used for microsatellite mining. SSRs were considered to contain motifs of two to six nucleotides and a minimum of five contiguous repeat units. Based on the MISA results, Primer 6.0 was used with the default settings to design primer pairs that would generate PCR products ranging from 100 to 280 bp in length. A total of 100 pairs of primers were designed (Additional file 1: Table S1) and validated by PCR in six P. canaliculata samples, including Guangzhou 1–2 (GZ1-2), Hainan1-2 (HN1-2), and Shaoguang1-2 (SG1-2) that were collected from three major invasive regions in Guangdong Province, China. PCR amplification was carried out as follows: an initial denaturation at 94°C for 4 min; 33 cycles of 94°C for 1 min (denaturation), 50°C for 30 s (annealing), and 72°C for 45 s (extension); followed by a final extension at 72°C for 8 min. The PCR products were analyzed by electrophoresis on a 8.0% non-denaturing polyacrylamide gel and silver stained.