A healthy individual of C. chinensis (weight 20.30 g, shell height 46.05 mm, shell width 33.27 mm) was collected from the rice snail breeding demonstration base in Ligao Village (23.37°N, 111.29°E), Liuzhou, Guangxi, China. Eight tissues were collected, namely the foot, intestine, hepatopancreas, gill, gonad, esophagus, blood, and lip. The samples were frozen with liquid nitrogen and stored at −80°C for later use. The trial was performed in accordance with the Institutional Animal Care and Use Committees of Guangxi Academy (CGA-00927).

The total RNA from these tissues was extracted separately, and equal amounts of RNA were mixed to construct a library after quality testing. Briefly, total RNA was extracted by grinding each tissue in TRIzol reagent (Life technologies, United States) on dry ice and processed following the manufacturer’s protocol. The integrity of the RNA was determined with an Agilent 2,100 Bioanalyzer and agarose gel electrophoresis. The purity and concentration of the RNA were determined with a Nanodrop micro-spectrophotometer (Thermo Fisher).

Referring to the methods of Huang et al., (2021). Briefly, mRNA was enriched using Oligo (dT) magnetic beads (PacBio biosciences, United States). The enriched mRNA (2 μg) was reverse transcribed into cDNA using a Clontech SMARTer PCR cDNA Synthesis Kit (Clontech, United States). PCR cycle optimization was used to determine the optimal amplification cycle number for the downstream large-scale PCR reactions. The optimized cycle number was then used to generate double-stranded cDNA. In addition, >5-kb size selection was performed using the BluePippin™ Size-Selection System and the resulting cDNA was mixed equally with non-size-selected cDNA. Large-scale PCR was performed for subsequent SMRTbell library construction. The cDNA samples were barcoded, then pooled and constructed as a “single” sample in the SMRTbell library, with different barcodes distinguishing the data from the different samples. The cDNAs were DNA damage repaired, end repaired, and ligated to sequencing adapters. The SMRTbell template was annealed to the sequencing primer and bound to polymerase, and sequenced on the PacBio Sequel II platform by Gene Denovo Biotechnology Co. (Guangzhou, China). A total of 57 Gb reads were generated from the SMRT library.

The raw sequencing reads of the cDNA libraries were analyzed using the unigene sequencing (Iso-Seq) pipeline supported by Pacific Biosciences (Gordon et al., 2015). First, high-quality circular consensus sequences (CCSs) were extracted from the subreads BAM file. The integrity of the transcripts was judged on the basis of whether the CCS reads contained 5′ primer, 3′ primer, and poly A structures. The sequences containing all three structures were called full-length sequences (FL reads). Subsequently, primers, barcodes, poly A tails, and concatemers of full passes were removed to obtain full-length non-chimeric (FLNC) reads. The FLNC reads were clustered to generate complete unigenes. Similar FLNC reads were clustered hierarchically using minimap2 to obtain a consistent sequence (unpolished consensus unigenes) (Li., 2018). The quiver algorithm was then used to further correct the consistent sequences. Two strategies were used to improve the accuracy of PacBio reads. First, we used non-full-length reads to polish the obtained cluster consensus Unigenes with the Quiver software to produce FL, polished, high-quality consensus sequences (accuracy of ≥99%). Second, Illumina short reads were obtained from the same samples with the LoRDEC to further correct low-quality isoforms (Salmela & Rival, 2014). On the basis of the results, both high-quality unigenes (prediction accuracy was ≥ 0.99) were combined for subsequent analysis.

To annotate them, the unigenes were BLAST analyzed against the NCBI non-redundant protein (Nr) database (https://www.ncbi.nlm.nih.gov), the Swiss-Prot protein database (https://www.expasy.ch/sprot), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa & Goto, 2000), and the COG/KOG database (https://www.ncbi.nlm.nih.gov/COG) with the BLASTx program (https://www.ncbi.nlm.nih.gov/BLAST/) at an E-value threshold of 1e-5 to evaluate sequence similarity to genes of other species. Gene Ontology (GO) annotation was performed using the Blast2GO software with the Nr annotation results for the unigenes (Conesa et al., 2005). Unigenes with the top 20 highest scores and no shorter than 33 HSP (High-scoring Segment Pair) hits were selected for Blast2GO analysis (Table 1). Then, functional classification of the unigenes was performed using the WEGO software (Ye et al., 2006).

Open reading frames (ORFs) were detected using the ANGEL software with the unigene sequences to obtain CDSs, protein sequences, and UTR sequences (Shimizu et al., 2006). Protein-coding sequences of the unigenes were aligned with hmmscan to Animal TFdb (https://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/) to predict TF families. CNCI (version 2) was used to assess the protein-coding potential of transcripts without annotations with default parameters for potential ncRNAs (Sun et al., 2013). ncRNAs were classified according to their secondary structures and sequence conservation.

The MIcroSAtellite (MISA, https://pgrc.ipk-gatersleben.de/misa/) tool was employed for microsatellite mining across the whole transcriptome. Initially SSRs of 1–6 nucleotide units were identified with the minimum repeat unit number defined as 15 for mono-nucleotides, 6 for di-nucleotides, 5 for tri-nucleotides, and 4 each for tetra-, penta-, and hexa-nucleotides. If the distance between two SSRs was shorter than 100 bp, they were considered one SSR.

To analyze AS events in the transcript unigenes, COding GENome reconstruction Tool (Cogent) was first used to partition transcripts into gene families based on k-mer similarity and reconstruct each family into a coding reference genome based on De Bruijn graph methods (Li et al., 2017). Then, the SUPPA tool was used to analyze AS events in the transcript unigenes (Alamancos et al., 2015).

Using the PacBio Sequel II sequencing platform, we obtained 57,428,331,313 bp of original transcript data (about 57 Gb) and a total of 26,658,437 subreads with an average subread length of 2,154 bp and an N50 length of 2,466 bp. After CCS analysis, 818,895 CCS reads were obtained with a mean CCS read length and pass number of 2,428 bp and 29, respectively. By clustering and correction, 697,861 FLNC reads were identified with an average length of 2,315 bp and an N50 length of 2,530 bp. Finally, the redundant reads were removed by CD-HIT and 26,312 unigenes with a mean length of 2,572 bp were obtained (Table 1). The length distribution of unigenes is shown in Figure 1. The majority of unigenes were between 1,000 and 5,000 bp, and the longest unigene was about 14,000 bp.

The 26,312 unigenes were annotated using four databases, with 24,151 (91.79%), 22,756 (86.49%), 13,934 (52.96%), and 16,423 (62.42%) unigenes matched to the Nr, KEGG, KOG, and Swiss-Prot databases, respectively (Figure 2A). A total of 24,169 (91,86%) unigenes were annotated in all the databases (Figure 2B). After aligning the unigene sequences with the Nr database using BlastX, the sequence with the best comparison result for each unigene was taken as the corresponding homologous sequence, the species to which the homologous sequence belonged was determined, and the homologous sequences of each species were counted. The results showed that the species with the highest homology was Pomacea canaliculata (17,089 unigenes), followed by Mizuhopecten yessoensis (1,199) and Aplysia californica (969) (Figure 2C).

A total of 13,214 sequences were annotated using the GO function classification library under the three broad categories: molecular function (17,422 unigenes), cellular component (48,494), and biological process (80,752). “Binding” (8,587, 49.29%), “cell” (8,587, 17.71%), and “cellular process” (10,215, 12.65%) were the most common terms among the annotated unigenes in the three categories mentioned above (Figure 2D).

In KEGG pathway analysis, the unigenes were assigned to five main categories: cellular processes (4,595 unigenes), environmental information processing (4,728), genetic information processing (1,797), metabolism (8,156), and organismal systems (10,363). The signal transduction (4,297, 16.33%) category had the largest number of unigenes. Next was infectious diseases (4,230, 16.08%) and cancers (3,068, 11.66%) (Figure 2E). Based on the KEGG annotation information we were able to further obtain Pathway annotations for unigene and focus on some of the more applicable pathways including PI3K-Akt signaling pathway (436, 5.34%), cGMP-PKG signaling pathway (389, 4.77%), Drug metabolism-cytochrome P450 (91, 1.12%), Platinum drug resistance (127, 1.56%), AMPK signaling pathway (267,3.27%), NF-kappa B signaling pathway (119, 1.46%), T cell receptor signaling pathway (71, 0.87%), Toll and Imd signaling pathway (80, 0.98%) (Table 2). From these, we also identified some candidate genes involved in drug metabolism and immune related, such as cytochrome P450 family (CYP1A1, CYP2J and CYP2U1), phosphoinositide-3-kinase regulatory subunit (PIK3), cGMP-inhibited 3′,5′- cyclic phosphodiesterase A (PDE3A), glutathione S-transferase (GST), 5′-AMP-activated protein kinase (PRKAG), NF-kappa-B inhibitor alpha (NFKBIA), cell division control protein 42 (CDC42), baculoviral IAP repeat-containing protein 2/3 (BIRC2_3), and others (Table 2). These provide data for the subsequent development of biologic drugs in snail.

All unigenes were classified into 25 categories in the KOG functional classification. (Figure 2F). The highest number of unigenes was annotated to the general function prediction only category (2,807, 10.67%), followed by the signal transduction mechanisms (2,563, 9.74%) and posttranslational modification, protein turnover, and chaperones (2,196, 8.35%) categories.