Sequencing of Transcriptome and Small RNA Revealed the Xenobiotic Metabolism-Related Genes and Potential Regulatory miRNA in Asian Tramp Snail

The Asian tramp snail, Bradybaena similaris (Ferusssac), is an invasive land snail species and has been a rising agricultural pest in south of China. As a pest, it also plays a role in transmission of Angiostrongylus cantonensis. However, present studies on this species are rare and the molecular information is limited. For this purpose, we sequenced the transcriptome and small RNA of B. similaris collected from citrus orchards. In total, 89,747 unigenes with an N50 size of 1287 bp and an average length of 817 bp were generated from ∼8.9 Gb transcriptome and 31 Mb clean reads were generated from ∼36 Mb small RNA library. To demonstrate the usefulness of these two datasets, we analyzed a series of genes associated with xenobiotic metabolism and core RNAi machinery. Analysis of the transcripts resulted in annotation of 126 putative genes encoding cytochrome P450 monooxygenases (CYP, 45), carboxyl/cholinesterases (CCE, 13), glutathione-S-transferases (GST, 24), and ATP-binding cassette transporters (ABC, 44). Analysis of the small RNA detected 42 miRNAs. In addition, four genes involved in small RNA pathways (miRNA, piRNA, and siRNA) were identified, and a total of 430 genes that can be targeted by miRNAs were predicted. Moreover, we found that a few miRNAs could target certain genes involved in xenobiotic metabolism. Therefore, we believe that these two datasets and the characterization of the identified/predicted genes will facilitate the molecular study of this species as well as other land snails with agricultural importance.

Introduction

The Asian tramp snail, Bradybaena similaris (Mollusca: Gastropoda: Pulmonata: Stylommatophora: Bradybaenidae), is an invasive land snail that has been a widely distributed pest. It usually inhabits garden plants, crops, and herbaceous grounds covered with organic debris. B. similaris is an increasingly important agricultural pest, feeding on various plants and being harmful to some economically important plants such as coffee and fruits, which has directly resulted in great economic losses (Serniotti et al., 2019). They start to come out looking for food at the beginning of March (temperature >10°C) and are most active from April to June (10–30°C) (Wang et al., 2016). In addition, this pest has also been reported as intermediate hosts of trematodes and nematodes (Alves et al., 2014; Morassutti et al., 2014), and it has been one of known intermediate hosts of Angiostrongylus cantonensis in China (Shan et al., 2008). Chemical molluscicides, such as niclosamide and metaldehyde, are usually applied to control snails (Lutinski et al., 2016; Wu et al., 2018). However, application of chemical agents to control pests is not sustainable due to low selectivity to target pest, resulting in pesticide resistance. Moreover, it can contaminate the environment and affect the beneficial animals (Lutinski et al., 2016; Tunholi et al., 2017). Recently, studies were reported showing that plant-derived molluscicides have obviously molluscicidal activity in the laboratory (Jia et al., 2019; Ke et al., 2019). It should be highlighted that the compounds present in the composition of plant-derived molluscicides are complex and some of them have high toxicity to non-target organisms, so evaluation and demonstration of the eco-toxicity information will be needed before using in the field (Pereira et al., 2020).

To better understand the molluscicidal activity of chemical pesticides or plant secondary metabolism, it is necessary to understand the molecular information of xenobiotic metabolism on snail. The molecular mechanism of xenobiotic detoxification can be divided into three phases (Van Leeuwen and Dermauw, 2016). In phase I, cytochrome P450 monooxygenases (CYPs) and carboxyl/cholinesterases (CCEs) can alter the structures of toxic compounds by introducing a nucleophilic functional group to produce a compound with higher activity and water solubility. Then, in phase II, the enzymes such as glutathione-S-transferases (GSTs) reinforce the water solubility by conjugation with endogenous molecules. In phase III, cellular transporters, such as ATP-binding cassette transporters (ABCs), transport all conjugates out of the cell. These genes involved in xenobiotic detoxification also showed a link between host plant adaptation and pesticides resistance (Dermauw et al., 2013). However, the molecular information of B. similaris is rare, which has hampered the identification of the genes involved in xenobiotic detoxification. In addition, these genes are also potential targets of RNA interference (RNAi) in enhancing toxicity of pesticides and host plant secondary metabolites (Niu et al., 2018).

RNA interference, based on post-transcriptional gene silencing, is a promising novel approach for controlling pest and a sustainable way in plant protection (Liu et al., 2020; Yan et al., 2020). However, the machinery of RNAi pathway in B. similaris is still unknown. RNAi is classified into three pathways, including small-interfering RNA (siRNA), microRNA (miRNA), and piwi-interacting RNA (piRNA). Briefly, for siRNA pathway, Dicer 2 cleavages the long dsRNA into 19- to 23-nt double-stranded siRNAs. These siRNAs then load to Argonaute 2 centered RNA-induced silencing complex (RISC), which finally inhibits the expression of the target gene. For miRNA pathway, primary RNAs are processed by Drosha and Pasha into shorter precursor miRNAs, then followed by the slicing of Dicer 1 into (21–24 nt) mature miRNAs. For piRNA pathway, the precursor piRNAs are cut into 26- to 32-nt piRNAs by Argonaute 3 to generate antisense piRNAs. These antisense piRNAs are trimmed by Aubergine and P-element induced wimpy testis (Piwi) and then generate sense piRNAs to create a feedback loop, known as ping-pong cycle (Dowling et al., 2016; Cooper et al., 2019).

In this study, we are aiming to provide the datasets of deep sequenced transcriptome and small RNA. The usefulness of these datasets will be demonstrated to identify genes involved in xenobiotic detoxification as well as the potential regulatory miRNAs. We believe that this study will facilitate the molecular study of this species as well as other land snails with agricultural importance.

Materials and Methods

Sample Collection

The Asian tramp snail B. similaris (Figure 1A) were collected from the citrus orchards in the Citrus Research Institute of Chinese Academy of Agricultural Sciences (CAAS), located on the Beibei District, Chongqing, China, in early spring (March) 2017. To avoid potential contamination from host plants, only the foot muscle from 10 adult snails were dissected for RNA extraction.

FIGURE 1

Figure 1. (A) Bradybaena similaris and citrus leaves. (B) Homology analysis of unigenes for B. similaris.

Transcriptome and Small RNA Sequencing

After washing with chilled PBS (pH 7.4) three times, the dissected foot muscles were homogenized in TRIzol reagent (Invitrogen, Carlsbad, CA, United States) to isolate RNA under the manufacturer’s protocol. About 1% agarose gel and a Nanodrop ONE spectrophotometer (GE Healthcare Bio-Science, Uppsala, Sweden) were used to confirm RNA quality and quantity, respectively. The library was constructed using a TruSeq Total RNA Sample Prep Kit (Illumina, San Diego, CA, United States). Then, sequencing was performed with PE150 in a HiSeq X Ten platform (Illumina). Low-quality reads and reads with adaptors and poly-N were removed from the raw reads. Then, these high-quality reads were analyzed with the Trinity program.¹ The obtained sequence processed by Trinity, which could not be extended on either end, was called unigene. After that, those assembled unigenes were annotated with a cutoff E-value of 10^–5 using BLASTx search against the following databases: NR; GO; COG/KOG/eggNOG; KEGG; Pfam; Swiss-Prot.

The small RNA library for B. similaris was constructed with small RNAs (18–30 nt) by PE50 in a HiSeq X Ten platform. Reads with adaptors, poly-N, low quality were removed. Then, clean reads were mapped to B. similaris transcriptome database using Bowtie software (version 1.0). Sequence alignment of the clean reads was carried out with the following database: Silva,² GtRNAdb,³ Rfam,⁴ and Repbase⁵ to filter ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and other ncRNA and repeats. Subsequently, the remaining clean reads were annotated by miRBase (V21) to detect known miRNAs. Novel miRNA secondary structures were predicted by Ranfold tools software (Version 2.0). MiRanda (Version 3.3) and RNAhybrid (Version 2.1) were used to predict miRNA targets.

Identification of Xenobiotic Metabolism Related Genes and Core RNAi Machinery

We retrieved transcripts encoding putative B. similaris genes related to xenobiotic metabolism (including CYPs, CCEs, GSTs, and ABCs) and core RNAi machinery (including miRNA, siRNA, and piRNA pathway) from B. similaris transcriptome. Firstly, BlastStation-Local 64 software was used to screen candidate unigenes associated with xenobiotic metabolism and candidate unigenes associated with core RNAi machinery. Known sequences from arthropods including Drosophila melanogaster, Tribolium castaneum, Apis mellifera, Bombyx mori, and Tetranychus urticae, and molluscs including Aplysia californica and Biomphalaria glabrata were used as queries for the homology searches. The matching sequences were then translated using the ExPASy-Translate tool⁶ and manually checked. The complete coding region was determined according to the open reading frame (ORF) finder.⁶ Finally, the selected candidate sequences of B. similaris with a minimum ORF length (CYPs: 200 nt; CCEs: 300 nt; GSTs: 200 nt) were aligned with known genes from other arthropods and molluscs using Mega 7.0 with the default settings in ClustalW for phylogenetic analysis (Kumar et al., 2016). For B. similaris ABCs, we used the nucleotide-binding domain (NBD) for phylogenetic analysis and extracted the NBDs using NCBI BLAST.⁷ A Neighbor-Joining analysis was performed using bootstrapping with 1000 replicates to construct the phylogenetic trees.

Results

Overview of Transcriptome

A total of 29,809,392 clean reads were generated (NCBI accession number SRR6981555) for B. similaris, encompassing about 8.91 Gb sequencing clean data with a GC content (guanine–cytosine content) of 44.97%. More than 89% clean reads showed quality scores higher than the Q30 (base calling accuracy of 99.9%). These clean reads were further assembled into 89,747 unigenes with an average length of 817 bp and an N50 size (at which 50% of assembly was covered) of 1287 bp (Table 1 and Supplementary Figure 1). Among 89,747 unigenes, a total of 26,331 unigenes were annotated against various databases (including NR, GO, COG, KOG, eggNOG, KEGG, Pfam, and Swiss-Prot) (Supplementary Table 1). Out of 89,747 unigenes, 23,481 unigenes (26.16%) were matched against the NCBI non-redundant (NR) protein database (Table 1 and Supplementary Table 1). According to the best hit in the NR protein database, 13,775 unigenes (58.73% of total NR annotation unigenes) had a strong homology with the California sea hare A. californica, followed by the owl limpet Lottia gigantea (1911 unigenes, 8.15%), the oyster Crassostrea gigas (948 unigenes, 4.04%), deposit feeder Capitella teleta (332 unigenes, 1.42%), the sea anemone Nematostella vectensis (275 unigenes, 1.17%), etc. (Figure 1B).

TABLE 1

Table 1. Summary of the transcriptome and small RNA libraries from Bradybaena similaris.

Based on the analysis of GO annotation, a total of 8240 unigenes were classified into cellular component, molecular function, and biological process (Supplementary Figure 2). The largest category was biological process (including 18 subcategories; 20,222 sequences, 48.48%), followed by the cellular component (including 18 subcategories; 11,644 sequences, 27.92%) and molecular function (including 16 subcategories; 9845 sequences, 23.60%) (Supplementary Table 2). A total of 7744 unigenes were grouped into 25 COG categories (Supplementary Figure 3). The “general function prediction” was the most common category (1563, 20.18%), followed by “translation, ribosomal structure and biogenesis” (791, 10.21%), and “posttranslational modification, protein turnover, chaperones” (674, 8.70%). A total of 10,908 B. similaris unigenes were mapped onto 273 KEGG pathways. The most abundant pathway was “ribosome” (469, 4.30%) from Genetic Information Processing, followed by “protein processing in endoplasmic reticulum” (289, 2.26%) from Genetic Information Processing and “purine metabolism” (247, 2.26%) from Metabolism (Supplementary Figure 4).

Identification of Genes Potentially Involved in Xenobiotic Metabolism and Core RNAi Machinery

Canonical xenobiotic metabolism genes including CYPs, CCEs, GSTs, and ABCs were analyzed based on the constructed transcriptome. In total, we identified 45 putative P450 genes. CYP2 clade has the majority P450 genes (29 genes), followed by CYP3 clade (eight genes), CYP4 clade (five genes), and Mitochondrial (three genes) (Figure 2 and Supplementary Table 3). A total of 13 putative CCEs genes were identified and grouped into four clades: J clade (four putative genes), N clade (one putative gene), M clade (one putative gene), and NC clade (seven putative genes) (Figure 3 and Supplementary Table 4). To better characterize GSTs in B. similaris, known genes from humans, plants, bacteria, and other animal species were also included to construct the phylogenetic tree. In total, 24 GST genes were assigned to seven classes: sigma (12 genes), pi (one gene), mu (three genes), zeta (two genes), omega (four genes), elongation factor 1-gamma (EF1Bgamma) (one gene), and theta (one gene) (Figure 4 and Supplementary Table 5). A total of 44 ABCs of B. similaris were assigned to subfamily A (nine genes), subfamily B (six genes), subfamily C (16 genes), subfamily E (three genes), subfamily F (five genes), and subfamily G (five genes) (Figure 5 and Supplementary Table 6).

FIGURE 2

Figure 2. Phylogenetic analysis of putative cytochrome P450 monooxygenase (CYP) genes. Neighbor-joining tree of B. similaris, Aplysia californica (Ac), Biomphalaria glabrata (Bg), Drosophila melanogaster (Dm), and Tetranychus urticae (Tu) CYP protein sequences. Numbers at each branch indicate the percentage of the times a node was supported in 1000 bootstrap pseudo-replications by the neighbor-joining method. The sequences used for phylogenetic tree construction are listed in Supplementary Table 3.

FIGURE 3

Figure 3. Phylogenetic analysis of putative carboxyl/cholinesterase (CCE) genes. Neighbor-joining tree of B. similaris, Aplysia californica (Ac), Biomphalaria glabrata (Bg), Drosophila melanogaster (Dm), and Bombyx mori (Bm) CCE protein sequences. Numbers at each branch indicate the percentage of the times a node was supported in 1000 bootstrap pseudo-replications by the neighbor-joining method. The sequences used for phylogenetic tree construction are listed in Supplementary Table 4.

FIGURE 4

Figure 4. Phylogenetic analysis of B. similaris glutathione S-transferases (GSTs). The phylogenetic tree was supported in 1000 bootstrap pseudo-replications by the neighbor-joining method. GenBank accession numbers for the sequences used in the phylogenetic tree are as follows: A. californica (Ac) Theta, XP_005089374; B. glabrata (Bg) Theta, XP_013081705; D. melanogaster (Dm)Theta, NP_610509; A. californica Mu, XP_005093728; B. glabrata Mu, XP_013097111; Thais clavigera Mu, ACD13785; T. urticae (Tu) Mu, XP_015781114; A. californica Pi, XP_005099402; B. glabrata Pi, XP_013095077; Homo sapiens (Hs) Omega, NP_004823; B. glabrata Omega, XP_013062909; D. melanogaster Omega, NP_648237; D. melanogaster Sigma, NP_523767; Bombyx mori Sigma, NP_001036994; D. melanogaster Delta, NP_524326; B. mori Delta, NP_001036974; D. melanogaster Epsilon, NP_611323; Anopheles gambiae (Ag) Epsilon, XP_319969; H. sapiens Alpha, NM145740; Rattus norvegicus Alpha, AAF37739; A. californica Alpha, XP_012945278; D. melanogaster Zeta, NP_649894; A. gambiae Zeta, XP_312009; Photobacterium ganghwense Beta, PSU09724; Rhizobiales bacterium Beta, WP_113235669; Arabidopsis thaliana Phi, NP_192161; Silene vulgaris Phi, AAA33931; A. thaliana Tau, NP_180506; Hordeum vulgare Tau, BAE44477; A. thaliana Lambda, NM_001125685; Pinus tabuliformis, Lambda, AGC13140; Pimephales promelas Rho, AAF78081; Oreochromis niloticus Rho, ABV24480; Glycine max EF1Bgamma, NP_001241223; B. glabrata EF1Bgamma, XP_013072347; H. sapiens EF1Bgamma, NP_001395; Mus musculus Kappa; and H. sapiens Kappa, AAS01180. The sequences used for phylogenetic tree construction are listed in Supplementary Table 5.

FIGURE 5

Figure 5. Phylogenetic analysis of B. similaris ABC protein NBDs. A set of ABC protein NBDs extracted from B. similaris, D. melanogaster (Dm), A. californica (Ac), and B. glabrata (Bg) were used into phylogenetic analysis. Numbers at each branch indicate the percentage of the times a node was supported in 1000 bootstrap pseudo-replications by the neighbor-joining method. The sequences used for phylogenetic tree construction are listed in Supplementary Table 6.

By using the homology genes in insects including T. castaneum, D. melanogaster, B. mori, and A. mellifera, we identified four core RNAi machinery, including two Aubergine from piRNA pathway, one Argonaute 1 from miRNA pathway, and one Dicer 2 from siRNA pathway from B. similaris (Supplementary Figure 5 and Table 7).

MiRNA Annotation and Potential Targets Prediction

The dataset generated 31,456,038 (86.56%) clean reads from 36,343,042 raw reads (NCBI accession number SRR6981556), of which 6,443,477 (17.73%) reads matched to reference transcriptome (Table 1). The sequence length of most small RNA was between 21 and 24 nt, of which sequence length in 22 nt was abundant (Figure 6A). MiRNAs were annotated based on the miRBase (version 21.0) database⁸; if the sequence completely matched a known miRNA, we then assigned the known miRNA name to the candidate B. similaris miRNA. If the sequence matched the seed regions (bases 2–8) but did not exactly match the remaining bases of the known miRNA, we then regarded it as novel miRNA. In total, we obtained 42 miRNA in B. similaris, including 19 known miRNAs (Table 2 and Figure 6B) and 23 novel miRNA (Supplementary Table 8 and Figure 6C). The preferred bases of B. similaris miRNA at the seed region positions that play an important role in miRNA complementary and target recognition were G(guanine), G, A(adenine)/U(uracil), U, U/G, A, and C(cytosine)/G (Figure 6D). The first base of the miRNAs (21–25 nt) indicated that U was the most used base while G was the least used base. The most common first base of the miRNAs with 22 nt was U (61.54%) (Figure 6E), which is the most preferred nucleotide for Dicer while recognizing and trimming the precursor miRNA.

FIGURE 6

Figure 6. Analysis of miRNAs length and base bias of B. similaris. (A) Length of all miRNAs. (B) Length of known miRNAs. (C) Length of novel miRNAs. (D) Analysis of the nucleotide bias at each position of the miRNAs. (E) Analysis of the nucleotide bias at the first position of the miRNAs.

TABLE 2

Table 2. List of known miRNAs from Bradybaena similaris.

To validate the possible biological context of the identified B. similaris miRNA for future study, we predicted the miRNA target from the small RNA dataset. A total of 91 genes and 341 genes could be targeted by 15 conserved known miRNA and 19 novel miRNA, respectively (Supplementary Tables 9, 10). Interestingly, certain gene xenobiotic metabolism can be targeted by a few novel miRNAs (Supplementary Table 11). One CYP gene (c131511.graph_c1) from CYP2 clade can be targeted by bsi-miR-novel 5, two ABCs (c122943.graph_c0, c98604.graph_c0) from ABC subfamily C can be targeted by bsi-miR-novel 6, and one CCE (c106892.graph_c0) from clade J can be targeted by bsi-miR-novel 8. This indicated that the miRNAs in B. similaris are potentially involved in the regulation of xenobiotic metabolism.

Discussion

Due to problems caused by intensive application of chemical pesticides, it is crucial to develop novel and effective alternative pest management strategies. RNAi is one of the most promising novel approaches for pest control, which triggers toxic effects in the target species without harming non-target organisms in the ecosystem (Christiaens et al., 2018). However, lack of genetic background for some pest, including the rising agricultural pest B. similaris in current research, has hampered in-depth studies and the development of better control strategies. Therefore, we investigated and described the results of a de novo deep sequencing of B. similaris transcriptome (8.91 Gb) and small RNAs (36.34 M) with Illumina sequencing technology.

In the transcriptome, a total of 89,747 unigenes were generated. The number of unigenes in our study is greater than those assembled in the marine gastropod Hemifusus tuba (61,575) (Li et al., 2019) but smaller than the snail Satsuma myomphala (103,744) (Kang et al., 2017), Echinolittorina malaccana (115,211) (Wang et al., 2014), and the rock scallop Crassadoma gigantea (96,006) (Cao et al., 2018). Mollusca is the second most species-rich phylum after Arthropoda in the animal kingdom, yet only a few genomic information has been published so far. For gastropods, the already sequenced genomes could be found in the sea hare A. californica (Moroz et al., 2006); the freshwater snail B. glabrata (Adema et al., 2017), Radix auricularia (Schell et al., 2017), and Lymnaea stagnalis (unpublished; genome data are available in the NCBI database: GCA_900036025.1); the sea snail Patella vulgata (Kenny et al., 2015), Conus tribblei (Barghi et al., 2016), and L. gigantea (Simakov et al., 2013); and the shelled pteropod Limacina bulimoides (Choo et al., 2020). In this study, out of 23,481 annotated unigenes against the NR protein database, 13,775 unigenes (58%) have strong homology with A. californica. Because these two species are distantly related organisms, the proportion is normal when compared with the sequenced transcriptome in other molluscs. The land snail S. myomphala showed 51% homologous unigenes with A. californica (Kang et al., 2017) and the freshwater mussel Cristaria plicata showed 40% homologous unigenes with the oyster C. gigas (Patnaik et al., 2016).

Xenobiotic metabolism genes on arthropods (i.e., insects and mites) have been well documented. However, they received little attention on molluscs and the majority of studies focused on the role in environmental pollutant metabolism, especially polycyclic aromatic hydrocarbons (Rewitz et al., 2006). For example, the snail Xeropicta derbentina was recommended as biomarkers monitoring pesticide exposure in the field because their B-type esterases are sensitive to organophosphorus insecticides (Laguerre et al., 2009). Based on the phylogenetic analysis, which is helpful for investigating the evolutionary relationships between genes among species, we identified a total of 126 genes related to potential xenobiotic metabolism, including 45 P450s, 13 CCEs, 24 GSTs, and 44 ABCs. According to the phylogenetic analysis, the total number of P450 unigenes in clade CYP2 (29 unigenes, 64%) is greater than other clades, which indicates that P450s in B. similaris were expanded in clade CYP2 when compared with B. glabrata (99 P450s, of which 33 in clade CYP2 and 35 in clade CYP3) (Adema et al., 2017) and arthropods, which have the majority of P450s in clade CYP3 and clade CYP4 such as mites (Bajda et al., 2015) and insects (Ai et al., 2011; Yu et al., 2015). Although most of them in B. similaris showed high homology with 2J2-like genes from A. californica and T. urticae, further studies are still needed to confirm this. Half of GSTs belong to Sigma subfamily, which may indicate that these genes in this subfamily are major detoxification genes. The clarification and characterization of these genes involved in the xenobiotic mechanism will be helpful to future studies in understanding the response of target genes to pesticides.

As a powerful gene function study tool, RNAi is widely used in insects that are also intensively explored to develop novel pest control strategy (Niu et al., 2018). Currently, RNAi are intensively investigated in invertebrates but relatively rare in molluscs (only 14 species reported) (Owens and Malham, 2015). Despite the limited studies, the successful case of RNAi-mediated gene silencing (by soaking) in the freshwater snail, B. glabrata, proved that RNAi is an achievable method for controlling snail pests (Knight et al., 2011). By using model insects as reference, we identified four genes associated with the core machinery RNAi pathways from the transcriptome of B. similaris, including one in miRNA pathway, one in siRNA pathway, and two in piRNA pathway. There were nine genes in B. glabrata, of which five were in miRNA pathway and four in piRNA pathway but no gene was reported in siRNA pathway (Adema et al., 2017). Although the full core RNAi machinery was not yet identified in B. similaris, it seemed that three RNAi pathways were conserved in this snail and provided confidence for further experimental assay to test the efficiency of RNAi in this species.

miRNA can regulate the post-transcription by repressing protein synthesis by either inhibiting the translation of mRNA or increasing mRNA degradation (O’Brien et al., 2018). We identified 42 miRNAs based on the constructed small RNA library, which is less than the number of miRNA in P. vulgata (45), L. gigantea (59), and B. glabrata (209) (Kenny et al., 2015; Adema et al., 2017). We found that a few xenobiotic mechanism genes can be targeted by miRNAs, which implied that miRNAs in B. similaris can regulate xenobiotic metabolism genes. Further experimental validations are required to confirm their involvement in pesticide detoxification.

Conclusion

In this study, we sequenced the transcriptome and small RNA of a rising agricultural pest B. similaris. These two datasets were useful in identifying molecular information for certain interesting genes, which were further demonstrated for the identification of the genes potentially involved in the xenobiotic metabolism and core RNAi machinery. In addition, miRNAs of this snail were also investigated, and the prediction of miRNAs in the regulation of xenobiotic metabolism genes was analyzed. However, these two datasets only serve as the starting point of molecular research in B. similaris; further characterization of xenobiotic metabolism genes is still required by experimental approach, such as RNAi-based gene function study. We believe that these two datasets could not only facilitate further research of the molecular understanding of B. similaris but also be used as a resource for the development of novel snail pest control approach.

Data Availability Statement

The datasets generated for this study can be found in the SRA (SRR6981555 and SRR6981556).

Author Contributions

QY and JN conducted the design. QY and BD contributed to the experiments. QY, WY, and FS conducted the statistical analysis. QY drafted the manuscript. JN and JW contributed to proofreading and editing this written work. All authors approved the final draft of the manuscript.

Funding

This study was funded by the earmarked fund for Modern Agro-industry (Citrus) Technology Research System of China (CARS-26).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.595166/full#supplementary-material

Footnotes

1. ^ http://github.com/trinityrnaseq/trinitynaseq/wiki
2. ^ http://www.arb-silva.de/
3. ^ http://lowelab.ucsc.edu/GtRNAdb/
4. ^ http://rfam.xfam.org/
5. ^ http://www.girinst.org/repbase/
6. ^ https://web.expasy.org/translate/
7. ^ https://blast.ncbi.nlm.nih.gov/Blast.cgi
8. ^ http://www.mirbase.org/

References

Adema, C. M., Hillier, L. W., Jones, C. S., Loker, E. S., Knight, M., Minx, P., et al. (2017). Whole genome analysis of a schistosomiasis-transmitting freshwater snail. Nat. Commun. 8:15451. doi: 10.1038/ncomms15451

PubMed Abstract | CrossRef Full Text | Google Scholar

Ai, J., Zhu, Y., Duan, J., Yu, Q., Zhang, G., Wan, F., et al. (2011). Genome-wide analysis of cytochrome P450 monooxygenase genes in the silkworm, Bombyx mori. Gene 480, 42–50. doi: 10.1016/j.gene.2011.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Alves, H., Tunholi-Alves, V. M., Tunholi, V. M., Golo, P., Bittencourt, V. R., and Pinheiro, J. (2014). Changes in the lipid profile of Bradybaena similaris (Ferussac, 1821) (Gastropoda, Xanthonychidae) during the development of Eurytrema coelomaticum (Giard and Billet, 1892) (Digenea, Dicrocoeliidae). Exp. Parasitol. 144, 52–56. doi: 10.1016/j.exppara.2014.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Bajda, S., Dermauw, W., Greenhalgh, R., Nauen, R., Tirry, L., Clark, R. M., et al. (2015). Transcriptome profiling of a spirodiclofen susceptible and resistant strain of the European red mite Panonychus ulmi using strand-specific RNA-seq. BMC Genom. 16:974. doi: 10.1186/s12864-015-2157-2151

CrossRef Full Text | Google Scholar

Barghi, N., Concepcion, G. P., Olivera, B. M., and Lluisma, A. O. (2016). Structural features of conopeptide genes inferred from partial sequences of the Conus tribblei genome. Mol. Genet. Genomics 291, 411–422. doi: 10.1007/s00438-015-1119-1112

CrossRef Full Text | Google Scholar

Cao, S., Zhu, L., Nie, H., Yin, M., Liu, G., and Yan, X. (2018). De novo assembly, gene annotation, and marker development using Illumina paired-end transcriptome sequencing in the Crassadoma gigantea. Gene 658, 54–62. doi: 10.1016/j.gene.2018.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Choo, L. Q., Bal, T. M. P., Choquet, M., Smolina, I., Ramos-Silva, P., Marletaz, F., et al. (2020). Novel genomic resources for shelled pteropods: a draft genome and target capture probes for Limacina bulimoides, tested for cross-species relevance. BMC Genomics 21:11. doi: 10.1186/s12864-019-6372-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Christiaens, O., Tardajos, M. G., Martinez Reyna, Z. L., Dash, M., Dubruel, P., and Smagghe, G. (2018). Increased RNAi efficacy in Spodoptera exigua via the formulation of dsRNA with guanylated polymers. Front. Physiol. 9:316. doi: 10.3389/fphys.2018.00316

PubMed Abstract | CrossRef Full Text | Google Scholar

Cooper, A. M., Silver, K., Zhang, J., Park, Y., and Zhu, K. Y. (2019). Molecular mechanisms influencing efficiency of RNA interference in insects. Pest Manag. Sci. 75, 18–28. doi: 10.1002/ps.5126

PubMed Abstract | CrossRef Full Text | Google Scholar

Dermauw, W., Wybouw, N., Rombauts, S., Menten, B., Vontas, J., Grbic, M., et al. (2013). A link between host plant adaptation and pesticide resistance in the polyphagous spider mite Tetranychus urticae. Proc. Natl. Acad. Sci. U.S.A. 110, E113–E122. doi: 10.1073/pnas.1213214110

PubMed Abstract | CrossRef Full Text | Google Scholar

Dowling, D., Pauli, T., Donath, A., Meusemann, K., Podsiadlowski, L., Petersen, M., et al. (2016). phylogenetic origin and diversification of RNAi pathway genes in insects. Genome Biol. Evol. 8, 3784–3793. doi: 10.1093/gbe/evw281

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, T. W., Wang, W., Sun, L. P., Lv, S., Yang, K., Zhang, N. M., et al. (2019). Molluscicidal effectiveness of Luo-Wei, a novel plant-derived molluscicide, against Oncomelania hupensis, Biomphalaria alexandrina and Bulinus truncatus. Infect. Dis. Poverty 8:27. doi: 10.1186/s40249-019-0535-537

CrossRef Full Text | Google Scholar

Kang, S. W., Patnaik, B. B., Hwang, H. J., Park, S. Y., Chung, J. M., Song, D. K., et al. (2017). Sequencing and de novo assembly of visceral mass transcriptome of the critically endangered land snail Satsuma myomphala: annotation and SSR discovery. Comp. Biochem. Phys. D 21, 77–89. doi: 10.1016/j.cbd.2016.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ke, W., Tu, C., Cao, D., Lin, X., Sun, Q., and Zhang, Q. (2019). Molluscicidal activity and physiological toxicity of quaternary benzo[c]phenanthridine alkaloids (QBAs) from Macleaya cordata fruits on Oncomelania hupensis. PLoS Neglect. Trop. Dis. 13:e0007740. doi: 10.1371/journal.pntd.0007740

PubMed Abstract | CrossRef Full Text | Google Scholar

Kenny, N. J., Namigai, E. K., Marletaz, F., Hui, J. H., and Shimeld, S. M. (2015). Draft genome assemblies and predicted microRNA complements of the intertidal lophotrochozoans Patella vulgata (Mollusca, Patellogastropoda) and Spirobranchus (Pomatoceros) lamarcki (Annelida, Serpulida). Mar. Genomics 24, 139–146. doi: 10.1016/j.margen.2015.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Knight, M., Miller, A., Liu, Y., Scaria, P., Woodle, M., and Ittiprasert, W. (2011). Polyethyleneimine (PEI) mediated siRNA gene silencing in the Schistosoma mansoni snail host, Biomphalaria glabrata. PLoS Neglect. Trop. Dis. 5:e1212. doi: 10.1371/journal.pntd.0001212

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S., Stecher, G., and Tamura, K. (2016). MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874. doi: 10.1093/molbev/msw054

PubMed Abstract | CrossRef Full Text | Google Scholar

Laguerre, C., Sanchez-Hernandez, J. C., Kohler, H. R., Triebskorn, R., Capowiez, Y., Rault, M., et al. (2009). B-type esterases in the snail Xeropicta derbentina: an enzymological analysis to evaluate their use as biomarkers of pesticide exposure. Environ. Pollut. 157, 199–207. doi: 10.1016/j.envpol.2008.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Bekaert, M., Wu, L., Mu, C., Song, W., Migaud, H., et al. (2019). Transcriptomic analysis of marine gastropod hemifusus tuba provides novel insights into conotoxin genes. Mar. Drugs 17:md17080466. doi: 10.3390/md17080466

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Jaouannet, M., Dempsey, D. A., Imani, J., Coustau, C., and Kogel, K. H. (2020). RNA-based technologies for insect control in plant production. Biotechnol. Adv. 39:107463. doi: 10.1016/j.biotechadv.2019.107463

PubMed Abstract | CrossRef Full Text | Google Scholar

Lutinski, J. A., Bonfante, F. P., Busato, M. A., and Guarda, C. (2016). Plant extracts, metaldehyde and saline solutions on the population control of Bradybaena similaris. Ciênc. Rural 46, 1316–1321. doi: 10.1590/0103-8478cr20151140

CrossRef Full Text | Google Scholar

Morassutti, A. L., Thiengo, S. C., Fernandez, M., Sawanyawisuth, K., and Graeff-Teixeira, C. (2014). Eosinophilic meningitis caused by Angiostrongylus cantonensis: an emergent disease in Brazil. Mem. Inst. Oswaldo Cruz 109, 399–407. doi: 10.1590/0074-0276140023

PubMed Abstract | CrossRef Full Text | Google Scholar

Moroz, L. L., Edwards, J. R., Puthanveettil, S. V., Kohn, A. B., Ha, T., Heyland, A., et al. (2006). Neuronal transcriptome of Aplysia: neuronal compartments and circuitry. Cell 127, 1453–1467. doi: 10.1016/j.cell.2006.09.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, J., Taning, C. N. T., Christiaens, O., Smagghe, G., and Wang, J.-J. (2018). Rethink RNAi in insect pest control: challenges and perspectives. Adv. Insect Physiol. 55, 1–17. doi: 10.1016/bs.aiip.2018.07.003

CrossRef Full Text | Google Scholar

O’Brien, J., Hayder, H., Zayed, Y., and Peng, C. (2018). Overview of MicroRNA biogenesis, mechanisms of actions, and circulation. Front. Endocrinol. 9:402. doi: 10.3389/fendo.2018.00402

PubMed Abstract | CrossRef Full Text | Google Scholar

Owens, L., and Malham, S. (2015). Review of the RNA interference pathway in molluscs including some possibilities for use in bivalves in aquaculture. J. Mar. Sci. Eng. 3, 87–99. doi: 10.3390/jmse3010087

CrossRef Full Text | Google Scholar

Patnaik, B. B., Wang, T. H., Kang, S. W., Hwang, H. J., Park, S. Y., Park, E. B., et al. (2016). Sequencing, De novo assembly, and annotation of the transcriptome of the endangered freshwater pearl bivalve, Cristaria plicata, provides novel insights into functional genes and marker discovery. PLoS One 11:e0148622. doi: 10.1371/journal.pone.0148622

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereira, L., Ribeiro, E. C. G., Brito, M. C. A., Silveira, D. P. B., Araruna, F. O. S., Araruna, F. B., et al. (2020). Essential oils as molluscicidal agents against schistosomiasis transmitting snails - a review. Acta Trop. 209:105489. doi: 10.1016/j.actatropica.2020.105489

PubMed Abstract | CrossRef Full Text | Google Scholar

Rewitz, K. F., Styrishave, B., Lobner-Olsen, A., and Andersen, O. (2006). Marine invertebrate cytochrome P450: emerging insights from vertebrate and insects analogies. Comp. Biochem. Physiol. C Toxicol. Pharmacol. 143, 363–381. doi: 10.1016/j.cbpc.2006.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Schell, T., Feldmeyer, B., Schmidt, H., Greshake, B., Tills, O., Truebano, M., et al. (2017). An annotated draft genome for Radix auricularia (Gastropoda, Mollusca). Genome Biol. Evol. 9:evx032. doi: 10.1093/gbe/evx032

PubMed Abstract | CrossRef Full Text | Google Scholar

Serniotti, E. N., Guzmán, L. B., Beltramino, A. A., Vogler, R. E., Rumi, A., and Peso, J. G. (2019). New distributional records of the exotic land snail Bradybaena similaris (Férussac, 1822) (Gastropoda, Bradybaenidae) in Argentina. Bioinvas. Rec. 8, 301–313. doi: 10.3391/bir.2019.8.2.12

CrossRef Full Text | Google Scholar

Shan, L., Yi, Z., Peter, S., and Xiao Nong, Z. (2008). Emerging angiostrongyliasis in mainland China. Emerg. Infect. Dis. 14, 161–164.

Google Scholar

Simakov, O., Marletaz, F., Cho, S. J., Edsinger-Gonzales, E., Havlak, P., Hellsten, U., et al. (2013). Insights into bilaterian evolution from three spiralian genomes. Nature 493, 526–531. doi: 10.1038/nature11696

PubMed Abstract | CrossRef Full Text | Google Scholar

Tunholi, V. M., Tunholi Alves, V. M., Monteiro, C. O., Silva, L. C. D., Dolinski, C. M., Castro, R. N., et al. (2017). Biological, biochemical and histological features of Bradybaena similaris (Gastropoda: Pulmonata) infected by Heterorabditis indica (Rhabditida: Heterorhabditidae) strain LPP1. Exp. Parasitol. 179, 28–35. doi: 10.1016/j.exppara.2017.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Leeuwen, T., and Dermauw, W. (2016). The molecular evolution of xenobiotic metabolism and resistance in chelicerate mites. Annu. Rev. Entomol. 61, 475–498. doi: 10.1146/annurev-ento-010715-023907

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Hui, J. H., Chan, T. F., and Chu, K. H. (2014). De novo transcriptome sequencing of the snail Echinolittorina malaccana: identification of genes responsive to thermal stress and development of genetic markers for population studies. Mar. Biotechnol. 16, 547–559. doi: 10.1007/s10126-014-9573-9570

CrossRef Full Text | Google Scholar

Wang, X., Cui, L., Zhai, J., and Ding, Z. (2016). Stable and clumped isotopes in shell carbonates of land snails Cathaica sp. and Bradybaena sp. in north China and implications for ecophysiological characteristics and paleoclimate studies. Geochem. Geophys. Geosyst. 17, 219–231. doi: 10.1002/2015gc006182

CrossRef Full Text | Google Scholar

Wu, C. H., Ma, Y., Zhang, Z. Q., Tang, K., Du, G. L., Wang, J., et al. (2018). Effect and cost-effectiveness of three commonly used molluscicides in largescale field application. Chin. J. Schistosom. Control 30, 619–624.

Google Scholar

Yan, S., Ren, B., Zeng, B., and Shen, J. (2020). Improving RNAi efficiency for pest control in crop species. Biotechniques 68, 283–290. doi: 10.2144/btn-2019-2171

CrossRef Full Text | Google Scholar

Yu, L., Tang, W., He, W., Ma, X., Vasseur, L., Baxter, S. W., et al. (2015). Characterization and expression of the cytochrome P450 gene family in diamondback moth, Plutella xylostella (L.). Sci. Rep. 5:8952. doi: 10.1038/srep08952

PubMed Abstract | CrossRef Full Text | Google Scholar