Transcriptional landscape of small non-coding RNAs reveals diversity of categories and functions in molluscs

ABSTRACT

Small non-coding RNAs (sncRNAs) are non-coding RNA molecules that play various roles in metazoans. Among the sncRNAs, microRNAs (miRNAs) guide post-translational gene regulation during cellular development, proliferation, apoptosis, and differentiation, while PIWI-interacting RNAs (piRNAs) suppress transposon activity to safeguard the genome from detrimental insertion mutagenesis. While an increasing number of piRNAs are being identified in the soma and germlines of various organisms, they are scarcely reported in molluscs. To unravel the small RNA (sRNA) expression patterns and genomic function in molluscs, we generated a comprehensive sRNA dataset by sRNA sequencing (sRNA-seq) of eight mollusc species. Abundant miRNAs were identified and characterized in all investigated molluscs, and ubiquitous piRNAs were discovered in both somatic and gonadal tissues in six of the investigated molluscs, which are more closely associated with transposon silencing. Tens of piRNA clusters were also identified based on the genomic mapping results, which varied among different tissues and species. Our dataset serves as important reference data for future genomic and genetic studies on sRNAs in these molluscs and related species, especially in elucidating the ancestral state of piRNAs in bilaterians.

KEYWORDS:

• miRNAs
• piRNAs
• gene regulation
• transposon silencing
• molluscs

Introduction

Small non-coding RNAs (sncRNAs) (<200 nucleotides (nt) long) are usually non-coding RNA molecules that play a crucial role in gene regulation at the transcriptional and post-transcriptional levels. They are mainly classified as either microRNAs (miRNAs), PIWI-interacting RNAs (piRNAs), or small interfering RNAs (siRNAs) [1,2]. Among these, only miRNAs have been extensively investigated in animals, while both piRNAs and siRNAs have not yet been sufficiently explored. The best-known function of miRNAs is regulation of gene expression, which affects development, proliferation, differentiation, apoptosis, oncogenesis, and tumorigenesis [3]. In molluscs, miRNAs have been identified and characterized in Pinctada fucata [4], Crassostrea gigas [5], Chlamys farreri [6], Littorina littorea [7], Mytilus galloprovincialis [8], and Hyriopsis cumingii [9]. To date, only five miRNAs from Haliotis rufescens and 64 miRNAs from Lottia gigantea have been deposited in miRBase (Version 22.0) [10]. A genome-wide profiling initiative has also been undertaken to discover miRNAs in molluscs [11]. However, the research on the identification and characterization of miRNAs in molluscs is far from complete.

piRNAs are a widespread mechanism to effectively suppress transposable element (TE) activity to safeguard the genome from detrimental insertion mutagenesis [12], a common issue in most animal germ lines. Interestingly, recent findings further imply that PIWI and piRNAs play vital roles in somatic cells [13]. Numerous piRNAs have been identified and characterized from both the soma and germlines of Nematoda [14], Porifera [15], Platyhelminthes [16], Annelida [17], Cnidaria [18,19], Echinodermata [17], Arthropods [20], and Chordata [21–23], whereas only a few piRNAs have been detected in molluscs [24–27]. In fact, considering the somatic expression of piRNAs in many bilaterians, it is parsimonious to assume that somatic piRNA expression represents the ancestral state that was established in an early protostomian ancestor [20,24]. The development of next-generation sequencing (NGS) has made it possible to identify high-throughput miRNAs and piRNAs in non-model organisms by small RNA sequencing (sRNA-seq).

Mollusca is the second largest phylum within Metazoa and is widely distributed throughout limnetic, oceanic, and intertidal aquatic ecosystems. Molluscs are remarkably resilient to harsh environmental conditions, such as wide fluctuations in salinity, temperature, and hydraulic pressure, as well as microbe-rich environments and prolonged exposure to air [11]. Recent advances in genomic studies of molluscs have provided insights into their biomineralization and immune and stress responses [28–30]. Although thousands of genes and numerous genomes have been discovered in molluscs [11], only a limited number of studies have addressed the diversity and function of sncRNAs. piRNAs were first identified in the somatic and gonadal tissues of the molluscs pacific oyster C. gigas and the great pond snail Lymnaea stagnalis [24], and subsequently detected in the pearl oyster P. fucata [25] and Biomphalaria glabrata snail [26]. As limited somatic piRNAs have been detected in molluscs thus far, it is difficult to conclude whether the observed level of somatic piRNA expression represents the ancestral state.

In the present study, we constructed a total of 60 sRNA sequencing libraries from somatic and gonadal tissues of eight molluscs to identify novel miRNAs and putative piRNAs. Tens of known and novel miRNAs were identified and annotated in each species, enriching the research on miRNAs in molluscs. Furthermore, abundant putative piRNAs were also discovered in the somatic and gonadal tissues of the rock oysters Crossostrea nippona, scallop Mizuhopecten yessoensis, red abalone Haliotis rufescens, orient clam Meretrix lamarckii, and mussel Mytilus galloprovincialis. The discovery of putative piRNAs in various molluscs may contribute to the elucidation of the ancestral state of piRNAs in bilaterians. This study presents new data about small RNAs that are useful for future miRNA and piRNA studies in molluscs and related species. These data would also serve as an important reference for genomic and genetic studies in molluscs.

Materials and methods

Sample collection

Eight species of molluscs (C. nippona, M. yessoensis, H. rufescens, Meretrix lamarckii, M. galloprovincialis, Sinotaia quadrata, Anadara broughtonii, and Neptunea arthritica) were collected from Tokyo Bay and the nearby fishery market in October 2020 and cultured separately within an aquarium for one week. Ten cultured specimens of each species were dissected, and tissues were collected from their adductor muscles (bivalves) or feet (gastropods), gills, mantles, and gonads. The harvested tissues were stored in 2 ml tubes containing 15× volumes of NucleoProtect® RNA solution (TaKaRa, Shiga, Japan) for 12 h at 4°C and then preserved at −80°C until further use.

RNA extraction, cDNA library construction and illumina sequencing

Total RNA was extracted from each somatic and gonadal tissue using the ReliaPrep^TM miRNA Cell and Tissue Miniprep System (Promega, WI, USA) according to the manufacturer’s instructions. The Qubit RNA Assay Kit in Qubit 3.0 Fluorometer (Life Technologies, CA, USA) and RNA ScreenTape Assay Kit in Agilent 2200 TapeStation (Agilent Technologies, Waldbronn, Germany) were used to assess the total RNA quantity and quality, respectively. Three isolated total RNA samples were mixed at equivalent concentrations to construct an sRNA-seq library for each tissue of the molluscs. Each tissue had two mixed libraries as replicates, except for the S. quadrata samples. sRNA-seq libraries were generated using the SMARTer smRNA-Seq Kit for Illumina (TaKaRa, CA, USA) following the manufacturer’s recommendations. The total RNA was first polyadenylated to provide a priming sequence for an oligo(dT) primer, and then the first strand of cDNA was synthesized by adding a second adapter sequence to the 3’ end. In the next step, full-length Illumina adapters were added during PCR amplification, including sequences required for clustering on an Illumina flow cell, Illumina TruSeq HT index sequences, and sequences bound by the Read Primer 1/2 sequencing primers. Library validation was performed using a High Sensitivity D1000 ScreenTape Kit in an Agilent 2200 TapeStation (Agilent Technologies). The cleanup products were purified, and size-selected using AMPure XP beads (Beckman Coulter, CA, USA) to obtain products <150 bp. The post-size-selection library was also validated using an Agilent 2200 TapeStation (Agilent Technologies), followed by sequencing on an Illumina HiSeq 2500 platform with a 50 bp single-end module (Macrogen, Tokyo).

Sequencing data pre-processing

To obtain high-quality sRNA reads, the raw data were filtered by removing: (1) reads containing a strong poly (A) tail, (2) reads containing poly (N), and (3) low-quality reads (average quality score ≤ 20). The 3 end adapter sequences were trimmed from sequencing reads using the fastx_toolkit (fastx_clipper – AAAAAAAA -c -l 15 -d 0 -Q 33) (http://hannonlab. cshl.edu/fastx_toolkit/index.html), followed by size filtration with a length range of 15–35 nt. The qualifying reads of these eight molluscs were then mapped to the reference genome sequences and Rfam database to remove the unmapped reads from the genome and rRNA, tRNA, snRNA, and snoRNA sequences using bowtie (v1.2.1, -f -k 3 -v2).

miRNA identification and annotation

The sRNA reads were pooled from different replicates and tissues for each species separately to perform miRNA prediction with miRBase 22.0 [10] and miRDeep2 (v.2.0.0.8) [31] using default parameters. The mapper.pl resource was used for preprocessing and genome mapping, miRDeep2.pl for de novo prediction, and quantifier.pl for expression profiling. Differentially expressed miRNAs were normalized to the total miRNA reads from each library between different tissues within the same species. The miRNA counts were normalized according to the reads per million (RPM) method.

miRNA target prediction

Based on the annotated genome results, 3’ UTR sequences were isolated from A. broughtonii and used for miRNA-target prediction. Target prediction for miRNAs was performed using miRanda [32] and RNA22 [33]. For the miRanda algorithm, an interacting score > 100 and a free energy <−10 kcal/mol were used to determine the sequence interactions between miRNAs and target genes. The RNA22 algorithm was used to determine the most favourable hybridization sites between miRNA and mRNA with a sensitivity of 63%, specificity of 61%, and seed size of 7, allowing a maximum of one unpaired base in the seed [11]. Genes that were predicted to be targeted by the same miRNAs using two different algorithms were considered to be target genes of specific miRNAs. Subsequently, the predicted target genes were assessed with gene ontology (GO) and KEGG databases using Blastx with an E-value <10⁻⁵. Blast2GO [34] was used to obtain GO annotation of the target genes based on Blastx hits. The network of miRNAs and their predicted target genes was drawn using NetworkAnalyst [35].

piRNA identification and processing

After miRNA identification, the miRNA reads were removed from the clean reads using Seqkit (v0.7.3, common module) [36]. Reads that did not match to any known ncRNAs were considered putative piRNAs and were subsequently identified as novel piRNAs. PPmeter [24], unitas [37], and proTRAC [38] were used to quantify and compare the amount of ongoing ping-pong amplification, search ping-pong signatures, calculate a Z-score for the enrichment of 10bp overlaps according to Zhang et al. [39], and predict genomic piRNA clusters using default parameters, respectively. Homologous piRNA clusters were detected using the NCBI BLASTn function. The relative expression patterns of piRNA clusters were calculated as mapped reads normalized for piRNA cluster length according to the transcripts per million reads (TPM) method [40].

Transposon element annotation

For TE annotation, we performed a de novo prediction of repetitive elements in the reference genome of each species using RepeatMasker (v4.0.7) [41] based on a de novo repeat library, constructed using RepeatModeler (v1.0.11) [42]. The analysis was conducted with the entire repeat dataset including the repeats localized in the predicted piRNA clusters of the eight molluscs. The enrichment of transposon categories was compared between the genome and piRNA clusters. A t-test was used to compare the transposon percentages among whole genome sequences and piRNA clusters. As no reference genome available, the genomes of C. gigas [43] were used as reference genome for C. nippona, Cyclina sinensis [44] for M. lamarckii, and Pomacea canaliculata [45] for S. quadrata, and N. arthritica, respectively.

Results

sRNA-seq in molluscs

A total of 262,882,185 raw reads were generated by Illumina sRNA-seq with a range of 21–46 million reads for each species (Table 1). After removing the low-quality and size-excluded reads, a total of 93,351,335 clean reads were obtained. In addition to the miRNA-sized (20–24 nt) reads detected, two different piRNA-sized reads (24–31 nt) were detected in the sRNA sequences of the molluscs. Two different size ranges of 25–27 nt and 29–31 nt were observed in the somatic and gonadal tissues of C. nippona, M. yessoensis, and M. lamarckii (Figure 1(a); Figure S1A,B), whereas dromedary piRNA-sized reads of 25–27 nt in length were detected in M. galloprovincialis, H. rufescens, S. quadrata, A. broughtonii, and N. arthritica (Figure 1(b), Figure S1C–F). The sequence composition of the clean reads displayed a first position nucleotide bias for uracil (U) (Figure 1(c)). sRNA identification and functional annotation results showed that non-coding RNAs, such as miRNA, snRNA, snoRNA, lncRNA, and coding nucleotides such as mRNA and rRNA were both detected by sRNA-seq (Figure 1(d)). However, a large proportion of the reads could not be annotated with the chosen reference of zebrafish (Danio rerio). In particular, numerous piRNA-sized reads could not be annotated since piRNA sequences are not conserved between species.

Figure 1. Characteristics of sRNAs in the investigated molluscs. (a) length distribution of clean reads in C. nippona; (b) length distribution of clean reads in M. galloprovincialis; (c) base composition at each site of the unique sRNA read from the 5’ terminal in C. nippona; (d) mapping statistics of clean reads among different RNA types in molluscs annotated by sRnatools using zebrafish (danio rerio) as the reference.

Table 1. Statistics of small RNA sequencing result of the eight investigated molluscs.

Species	Total Bases	Read Count	Clean Read	GC (%)	Q20%)	Q30%)
C. nippona	2,343,411,240	45,949,240	12,712,134	46.36	94.29	90.22
M. yessoensis	1,450,135,326	28,434,026	9,427,110	43.16	90.68	83.77
H. rufescens	1,678,385,214	32,909,514	11,882,093	44.10	92.21	87.04
M. lamarckii	1,511,284,632	29,633,032	10,477,600	45.31	95.05	91.44
M. galloprovincialis	1,531,547,289	30,030,339	11,546,665	43.44	93.15	88.54
S. quadrata	1,075,614,123	21,090,473	4,541,116	47.40	93.73	89.37
A. broughtonii	2,273,109,423	44,570,773	21,881,222	44.01	93.63	89.23
N. arthritica	1,543,504,188	30,264,788	10,883,395	44.86	93.07	88.34

miRNA identification in molluscs

The clean reads were used to identify miRNAs in molluscs. A total of 72, 85, 69, 42, 73, 57, 42, and 43 miRNAs were identified from the eight mollusc species (Table S1). Among these miRNAs, ten were detected in all eight molluscs (Figure 2(a)): miR-92a, miR-317, miR-7024-3p, miR-4524a, miR-96-5p, miR-87, miR-2424, miR-365, miR-8935 and miR-146b-3p. Another ten miRNAs: miR-1990c-3p, miR-7171-5p, miR-1985, miR-22-3p, miR-2a, miR-9007, miR-876, miR-12336-5p, miR-193a, and miR-283 were also detected in five or more of the mollusc species. Based on the miRNA expression patterns, principal component analysis (PCA) clustered C. nippona, M. lamarckii, S. quadrata, A. broughtonii, and N. arthritica together, but segregated M. yessoensis, H. rufescens, and M. galloprovincialis (Figure 2(b)). The expression patterns of miRNAs varied among different tissues, and differed significantly between somatic and gonadal tissues in A. broughtonii, H. rufescens, and S. quadrata (Figure 2(c); Figure S2), whereas the lengths and sequences of miRNAs were mostly conserved, especially in the seed region (2–8 bp). For instance, the mature sequence of miR-22-3p was conserved in all investigated mollusc species, with a very similar precursor sequence (Figure 2(d)). In addition, abundant isoforms were identified in numerous miRNAs in molluscs, as shown in Figure 2(e).

Figure 2. miRNA identification in the investigated molluscs. (a) venn diagram of common and unique miRnas identified in the investigated molluscs; (b) Principal component analysis (PCA) of miRNA expression patterns in molluscs; (c) boxplot of miRNA expression density in molluscs. Adductor muscle collected from bivalves and foot muscle from gastropods; (d) alignment of precursors of miR-22-3p in molluscs; (e) alignment of precursors of miR-87 in molluscs; i1/i2/i3 indicates multiple miRNA isoforms; ‘pre’ represents precursor. The underline shows mature miRNA sequences.

miRNA target prediction in molluscs

To further understand the function of miRNAs in molluscs, we performed miRNA target prediction using miRanda and RNA22 for A. broughtonii. More than ten thousand target sites were predicted to be targeted by A. broughtonii miRNAs (Figure 3(a)). The genes predicted to interact with the same miRNAs using different algorithms were considered as strong candidates for miRNA targets and selected for further functional analysis. A total of 1,918 target sites were anticipated to be miRNA-mRNA interactions among 56 miRNAs and 1,223 genes in A. broughtonii. Most miRNAs were predicted to regulate multiple genes, with an average of 34 target genes. miR-317 had the greatest number of predicted target genes, followed by miR-2a, miR-8935, miR-3885-3p and miR-4187-5p (Figure 3(b)). miRNAs form a network with target genes associated with a complex process (Figure 3(c)). The enrichment analysis showed that 120 GO terms (Table S2) and 38 KEGG pathways (Figure 3(d)) associated with metabolism, energy conversion, signal transduction, and multiple molecular and cellular processes were enriched in the predicted target genes of miRNA in A. broughtonii.

Figure 3. Target prediction of miRnas in A. broughtonii. (a) venn diagram of the number of predicted target sites between genes and miRNAs using miRanda and RNA22 algorithms. (b) the number of target genes for each miRNA. (c) miRnas-target genes regulatory network of miRnas in A. broughtonii. (d) kyoto encyclopaedia of genes and genomes (KEGG) pathways, terms, and statistics of the KEGG pathway enrichment analysis for the predicted target genes of miRnas in A. broughtonii.

Identification and characteristics of piRNAs in mollusks

After removing the known ncRNAs, including the identified miRNA reads, the remaining reads were considered putative reads for piRNA processing. Abundant piRNA-sized reads with lengths of 24–32 nt were observed in the somatic and gonadal tissues of the molluscs (Figure 4(a)). Thousands of piRNAs were identified from the cluster loci in the genome (Supplementary sequences). In C. nippona, a total of 46 piRNA clusters (271.08 kb, 0.055% of the reference genome of C. gigas, assigned with 16.29% piRNAs) were identified using the reference genome of C. gigas, 7 piRNA clusters (49.08 kb, 0.005% of the genome) in M. yessoensis, 14 piRNA clusters (113.15 kb, 0.008% of the genome) in H. rufescens, 1 piRNA cluster (7.35kb, 0.001% of the reference genome of C. sinensis, assigned with 1.19% piRNAs) in M. lamarckii, 5 piRNA clusters (36.77kb, 0.003% of the genome) in M. galloprovincialis, and 6 piRNA clusters (44.41kb, 0.005% of the genome) in A. broughtonii (Table 2). However, no homologous sequence was detected among the piRNA clusters identified in molluscs. In addition, some putative piRNA sequences were not found in the identified piRNA clusters in molluscs.

Figure 4. Characteristics of piRNAs in the investigated molluscs. (a) length distribution of filtered known ncRNAs; (b) ping-pong pair reads of putative piRNAs; (c) ping-pong Z₁₀ score for the enrichment of 10-bp overlaps; (d) ping-pong matrices illustrate frequent length-combinations of ping-pong pairs (sequences with a 10-bp 5’ overlap; E: putative piRNA length distribution and 1 U/10A bias for ping-pong reads.

Table 2. Identification of miRNAs and piRNAs in the investigated molluscs.

Species	Known miRNAs	Novel miRNAs	Unique piRNAs	piRNA clusters	Cluster length (kb)	Clustered piRNAs
C. nippona*	63	9	32,947	46	271.08	16.29%
M. yessoensis	70	15	9,114	7	49.08	4.92%
H. rufescens	63	6	17,637	14	113.15	6.64%
M. lamarckii*	37	5	896	1	7.35	1.19%
M. galloprovincialis	62	11	12,012	5	36.77	1.40%
A. broughtonii	49	8	10,159	6	44.41	11.83%
S. quadrata*	40	2	3	0	0	0
N. arthritica*	41	2	88	0	0	0

*Genomes of closely related species were used as references for piRNA identification.

A significant 10 nt overlap between sense and antisense strands was detected in most of the molluscs except S. quadrata and N. arthritica, due to the lack of a reference genome for read mapping (Figures 4(b,c)), which is consistent with piRNAs generated by the ping-pong amplification mechanism. In gonadal tissues, most ping-pong pairs predominantly combined piRNAs of length 29/30 nt (27/28 nt in H. rufescens), suggesting homotypic ping-pong amplification (Figure 4(d)). In somatic tissues, ping-pong pairs combined piRNAs 29/30 nt and 25/26 nt in length, suggesting both heterotypic and homotypic ping-pong amplification. Furthermore, the ping-pong pairs were heavily biased towards uridine (U) at the 5 position and adenine (A) at position 10 (Figure 4(e)).

All putative piRNA reads from each species were re-mapped with the identified piRNA clusters to calculate the relative expression performance using the TPM method. In C. nippona (Figure 5(a)), M. lamarckii (Figure 5(b)), and M. yessoensis (Figure 5(c)), the piRNA clusters were highly expressed in the gonad tissues, whereas piRNA clusters were highly expressed in the somatic tissues of A. broughtonii (Figure 5(d)), H. rufescens (Figure 5(e)), and M. galloprovincialis (Figure 5(f)). In addition, different individuals also showed different expression patterns in these piRNA clusters, as seen in the gill and foot muscle of individual 2 which showed high expression in H. rufescens, and the gill tissue of individual 1 and mantle tissues of individuals 1 and 2 which showed high expression in M. galloprovincialis.

Figure 5. Heat map of expression levels of piRNA clusters in molluscs somatic and gonadal tissues. (a) C. nippona; (b) M. lamarckii; (c) M. yessoensis; (d) A. broughtonii; (e) H. rufescens; (f) M. galloprovincialis. ad: adductor muscle; Fo: foot muscle; gi: gill tissue; ma: mantle tissue; go: gonad.

piRNA-mediated TE suppression

For TE annotation, the genomes of C. gigas and C. sinensis were used as reference genomes for C. nippona and M. lamarckii, respectively. The C. gigas genome contains 176.55 Mb repetitive sequences, accounting for 37% of the genome (Figure 6(a)); however, the identified piRNA clusters in C. nippona contained 22% transposon sequences (Figure 6(b)), which was lower than that of the reference genome of C. gigas. The piRNA clusters identified in M. lamarckii also contained a lower percentage of transposons compared to the reference genome of C. sinensis. However, the piRNA clusters identified in M. galloprovincialis, H. rufescens, M. galloprovincialis, and A. broughtonii had a higher percentage of transposons compared to reference genome sequences of the same species (Figures 6(d,e); Figure S3). Annotation with transposon sequences revealed that they were rich in most piRNA clusters compared to the whole genome landscape, except in C. nippona and M. lamarckii, using the closely related reference genomes of C. gigas and C. sinensis (Fig. S4). Interestingly, the composition of transposons in piRNA clusters does not reflect the transposon landscape of the whole genome. Instead, piRNA clusters were rich in LTR/Pao, LTR/Gypsy, and DNA/Zator in C. nippona (Figure 6(c)), and LINE/Penelope, DNA/Harbinger, LTR/Gypsy, RC/Helitron, LINE/L1, and SINE/tRNA in M. yessoensis (Figure 6(f)), and a dominant unknown transposon detected in other molluscs (Fig. S3).

Figure 6. Transposon annotation in genome and piRNA clusters in the investigated molluscs. (a) representation of transposons in the genome of C. gigas, plotted by divergence [%] from transposon consensus. (b) representation of transposons within piRNA clusters of C. nippona plotted by divergence [%] from transposon consensus. (c) prominent transposons that are enriched or depleted in C. nippona piRNA clusters. (d) representation of transposons in the genome of M. yessoensis, plotted by divergence [%] from transposon consensus. (e) representation of transposons within piRNA clusters of M. yessoensis, plotted by divergence [%] from transposon consensus. (f) prominent transposons that are enriched or depleted in M. yessoensis piRNA clusters.

Discussion

Mollusca is the second largest phylum after Arthropoda, however the study of small RNAs is far from complete in molluscs. In the present study, a total of 60 sRNA-seq libraries were constructed for small RNA discovery in molluscs. Tens of known and novel miRNAs were identified from each species, and the characteristics and expression patterns of these miRNAs were investigated in different tissues and species. The target genes of miRNAs were also predicted using data from A. broughtonii, which further improved the understanding of miRNAs in molluscs. In addition, abundant putative piRNAs were detected in both somatic and gonadal tissue in molluscs. The dataset obtained from the present study can serve as important reference data for future genomic and genetic studies on small RNAs in molluscs.

A genome-wide computational prediction of homologous miRNAs was performed in 35 species of molluscs [11], however it cannot reflect the situation of miRNA expression in vivo. In the present study, sRNA-seq was effectively used to detect ncRNAs in molluscs, enabling the discovery of known and novel miRNAs. With a limit of sequencing samples of tissue and development stage, it is impossible to discover every miRNA candidate, regardless of sequencing depth or which developmental stages are considered [46]. An alternative method is to complement miRNA identification with genome scanning to obtain more comprehensive and reliable miRNA expression profiles and functional analysis. It is also worth discovering miRNAs in molluscs from an evolutionary perspective, as molluscs are an evolutionarily important, yet relatively under-examined group.

Among the identified miRNAs, ten miRNAs were detected in all investigated molluscs, and ten miRNAs were found in at least five mollusc species. In fact, miR-2 is an invertebrate-specific miRNA family that has many family members and plays a vital role in the immune response [47]. miR-87 is a critical factor required in neurons to reactivate dendritic growth both in developmental remodelling and following injury in invertebrates [48]. In contrast, miR-1990c-3p was specifically discovered in molluscs [11], which were highly expressed in mantle tissues (Table S1). Matrix metalloproteinases (MMPs) were predicted to be targeted by miR-1990c-3p in P. fucata [4], which indicated that they might play a role in biomineralization by regulating the formation of matrix proteins. These highly conserved miRNAs have a long evolutionary history, suggesting that they may play crucial roles in essential biological processes in many organisms.

The prediction of miRNA targets is of central importance for the functional understanding of miRNAs. However, functional studies of specific miRNAs are still rare in molluscs. To date, miR-15b [49], miR-29a [50], miR-183 [51], and miR-2305 [52] have been shown to participate in nacre formation in molluscs, whereas miR-124 is associated with neurogenesis in molluscs [53,54]. In the present study, we performed an algorithmic prediction of miRNA targets in A. broughtonii using multiple bioinformatics tools. Thousands of targets were identified based on the physical sequence interaction between miRNA and mRNA 3’ UTR [55,56]. Computational prediction of miRNA targets always yields a large proportion of false positives, especially in unverified miRNAs [57]. Functional annotation and enrichment analyses of the predicted miRNA targets can provide a global understanding of miRNA function in a specific organism. A complex miRNA-mediated regulatory network was generated based on the predicted miRNA-mRNA interactions in A. broughtonii, but experimental work is required to confirm the authenticity of miRNA-targets. Hence, experimental verification is vital for target validation in this study.

In addition to the miRNAs detected in these molluscs, abundant piRNA-sized reads, 25–31 nt in length, were also discovered in the somatic and gonadal tissues of molluscs by sRNA-seq. To date, piRNAs have not yet been detected in Protozoa but have been in Nematoda [14,58], Porifera [15,17], and Chordata [21–23,59–61]. Intriguingly, ubiquitously expressed piRNAs were discovered in soma and germlines in most invertebrates, whereas very few piRNAs were discovered in the soma outside the germline in vertebrates. In molluscs, piRNAs were recently discovered to be ubiquitously expressed in both the somatic and gonadal tissues [24–27]. In the present study, we also performed piRNA discovery in diverse Mollusca species, which will provide fundamental knowledge to elucidate the ancestral state of piRNAs in bilaterians from an evolutionary perspective.

The piRNAs newly identified from molluscs exhibit characteristics such as uridine at the 5’ terminal (1 U) and a length of 24–32 nt, which is consistent with piRNAs in most animals [62]. Genomic mapping of piRNAs in Drosophila, zebrafish, mouse, and primates has shown that most piRNAs map within discrete genomic intervals, called piRNA clusters [59,63,64]. In the primary biogenesis pathway, long piRNA precursors are transcribed from specific genomic loci piRNA clusters, cleaved, and modified by complex factors in the cytoplasm to form the piRNA-induced silencing complex (piRISC) [65,66]. Subsequently, the functional piRISC binds to specific sequences to participate in transcriptional gene silencing of TE in the nucleus or post-transcriptional gene silencing of the transcript undergoes a ping-pong amplification loop in the cytoplasm [24,67,68]. A strong ping-pong amplification signal was detected through comprehensive computational analysis of the piRNA population of molluscs, which implies an important role of piRNA in post-transcriptional gene silencing in the cytoplasm.

In molluscs, tens of piRNA clusters with a size of 2.7 ~ 17 kb were identified based on genomic mapping data. The pronounced clustering suggests that piRNAs are processed from long primary transcripts, although the piRNA clusters constitute less than 0.1% of the genomes. In general, the piRNA clusters expressed were not substantially different among somatic tissues but varied significantly between somatic and gonadal tissues. The piRNAs and piRNA clusters also exhibited different expression patterns in the somatic and gonadal tissues of P. fucata [25,27] The tissular difference of piRNA expression might imply different functions of piRNAs in molluscs. The key function of piRNAs is to suppress transposon activity via transcriptional or post-transcriptional silencing mechanisms, thereby maintaining germline genomic integrity [69,70]. As TE annotation was limited in molluscs, we performed a de novo annotation of TE in the molluscs genomes however, a large number of identified transposons are not well classified. In addition, the reference genomes of molluscs have not been well annotated, regardless of the quality or quantity of the assembled genomes. To date, less than 50 mollusc genomes have been studied, which is a small proportion of total mollusc species [11]. Of the eight species in this study, four species lacked their own reference genome, therefore closely related genomes were used as references. However, piRNAs show high diversity among species and as such closely related reference genomes can be used for miRNA detection, but are not appropriate for novel piRNA analysis. The piRNA clusters identified from H. rufescens, M. galloprovincialis, and A. broughtonii were significantly enriched in TE compared to the reference genomes, indicating the key role of piRNAs in TE silencing in molluscs. However, a decreased proportion of TE was detected in the piRNA clusters of C. nippona and M. lamarckii compared with a closely related reference genome. In the future, the reference genomes of C. nippona and M. lamarckii will enhance the functional analysis of piRNAs in these two species.

In addition to transposon silencing, the PIWI-associated piRISCs also plays an important role in genome rearrangement [71], epigenetic regulation [72], gene regulation and development [73,74], and virus defence [75]. Metagenomic sequencing of small RNAs indicated the presence of endogenous RNA or DNA virus-derived piRNA expression in divergent animal phyla, including Cnidaria, Echinodermata, and Mollusca [17]. More evidence on endogenous viral element-derived piRNAs supports the hypothesis that they mediate antiviral immunity, similarly to CRISPR RNAs in prokaryotes [76,77]. It has also been found that somatic piRNAs can regulate endogenous gene expression via miRNA-like pairing rules to regulate gene expression, especially in arthropods [20,74]. A previous study also revealed a piRNA-mediated maternal mRNA decay during the maternal-to-zygotic transition in the mosquito Aedes [78]. Recent research also indicated that somatic piRNAs regulate endogenous gene expression in P. fucata somatic tissues [27]. piRNA-mediated gene regulation will provide comprehensive insights into post-transcriptional gene silencing in animals. The piRNAs might have a great diversity of functions in molluscs, including but not limited to transposon silencing. Further studies are needed to fully understand piRNAs in mollusc somatic and gonadal tissues.

In summary, the miRNAs, piRNAs, and the piRNA clusters defined in our study provide important insights into the evolutionary origins and functional understanding of small non-coding RNAs in molluscs. These findings are useful for the identification of miRNAs and piRNAs in molluscs and other closely related species and contribute to the expanding knowledge base for miRNA and piRNA analysis in molluscs as well as for elucidating the ancestral state of somatic piRNAs in bilaterians.

Abbreviations

miRNA: microRNA; piRNA: PIWI-interacting RNA; piRISC: piRNA-induced silencing complex; siRNA: small interfering RNA; sncRNA: small non-coding RNA; sRNA-seq: small RNA sequencing; TE: transposon element.

Author contributions

S.H. designed and performed the experiments, analysed the data and drafted the manuscript. K.Y. and S.K. collected the samples and assisted with the data analysis. S.A. conceived the study, supervised the work, and co-wrote the manuscript. All authors reviewed the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The raw sequencing reads from smRNA-seq are available in the DNA Data Bank of Japan (DDBJ) Sequence Read Archive (DRA) under the accession number DRA012198.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15476286.2024.2348893.

Additional information

Funding

This study was supported by grants from the National Key Research and Development Program of China [2018YFD0900601], the Japan Society for the Promotion of Science [JP24248034] and the Research Fellowship of the Japan Society for the Promotion of Science Postdoctoral Fellowship for Overseas Researchers [P20395].