A survey of miRNAs involved in biomineralization and shell repair in the freshwater gastropod Lymnaea stagnalis

MicroRNAs (miRNAs) are a deeply conserved class of small, single stranded RNA molecules that post-transcriptionally regulate mRNA levels via several targeted degradation pathways. They are involved in a wide variety of biological processes and have been used to infer the deep evolutionary relationships of major groups such as the Metazoa. Here we have surveyed several adult tissues of the freshwater pulmonate Lymnaea stagnalis (the Great Pond Snail) for miRNAs. In addition we perform a shell regeneration assay to identify miRNAs that may be involved in regulating mRNAs directly involved in the shell-forming process. From seven mature tissues we identify a total of 370 unique precursor miRNAs that give rise to 336 unique mature miRNAs. While the majority of these appear to be evolutionarily novel, most of the 70 most highly expressed (which account for 99.8% of all reads) share sequence similarity with a miRBase or mirGeneDB2.0 entry. We also identify 10 miRNAs that are differentially regulated in mantle tissue that is actively regenerating shell material, 5 of which appear to be evolutionarily novel and none of which share similarity with any miRNA previously reported to regulate biomineralization in molluscs. One significantly down-regulated miRNA is predicted to target Lst-Dermatopontin, a previously characterized shell matrix protein from another freshwater gastropod. This survey provides a foundation for future studies that would seek to characterize the functional role of these molecules in biomineralization or other processes of interest.


Introduction
The diversity of proteins that animals employ to biomineralize is steadily being revealed, and it is now clear that both deeply conserved and lineage specific elements play important roles in this process. Molluscs, a diverse clade of metazoans whose success is in part due to the evolutionary plasticity of the calcified shell, boast a diverse catalogue of shell-forming proteins that control all aspects of the biomineralization process from the initiation of mineral deposition, to regulation of the CaCO 3 polymorph and inhibition of crystal growth. While the hemolymph (the molluscan equivalent of blood) and neurons have been implicated in regulating shell formation [1][2][3] and patterning [4], the mantle is the primary tissue that co-ordinates shell formation. The mantle epithelium is a highly complex tissue that is innervated, muscularized and rich in secretory and sensory cells [5,6] and generally lies between the shell and the soft tissues of conchiferan (shelled) molluscs [7]. All molluscs secrete shell-forming proteins (and other biomolecules such as polysaccharides and lipids) from the mantle tissue. These biomolecules are secreted either into the mantle cavity (also known as the pallial space, terms used to refer to the volume between the mantle tissue and the shell) or as a water-insoluble complex (the periostracum) that forms a matrix upon and within which calcification proceeds [8].
While the inventories of molluscan shell-forming proteins continue to grow at an exponential rate [9,10], the ways in which these proteins are regulated and the underlying gene regulatory networks (GRNs) that have evolved to coordinate these processes [11] remain more obscure. Transcription factors and signaling molecules are arguably the primary elements of any GRN, and several have been directly implicated in coordinating the specification and differentiation of shell-forming cells in molluscs, and in directly regulating the biomineralization process [12][13][14][15][16][17]. Other regulatory molecules such as microRNAs (miRNAs) are also responsible for tuning mRNA transcripts to appropriate levels during many biological processes. The molecular machinery for the production of these small molecules is deeply conserved, and has been characterized for several molluscs, but with an emphasis on bivalves [18]. In addition, because the presence/absence of miRNAs has been used for deep phylogenetic reconstructions of animal relationships the miRNA complements of diverse metazoan genomes, including molluscs, have been surveyed [19][20][21][22][23].
A number of recent studies have identified miRNAs that play key roles in the process of biomineralization in molluscs [24][25][26][27][28][29][30][31], however all of these studies were focused on bivalves, in particular on the pearl forming genus Pinctada. In several cases these studies identify miRNAs that play deeply conserved roles in biomineralization. For example miR-29a and miR-183 participate in shell formation in P. martensii [26,27] and also play roles in vertebrate bone deposition [32][33][34]. Other studies have identified and classified the complete miRNA contents of gastropod genomes [35,36] while others have investigated the role of miRNAs in a variety of processes. For example, a conserved miRNA (miR-137) is associated with memory in the gastropod Lymnaea stagnalis [37], while miR-124 is associated with neuronal regeneration in the same species [38]. Biggar et al. [39] demonstrated an association between the upregulation of several conserved miRNAs and extended periods of freezing and anoxia in Littorina littorea. However no study has focused on the expression of miRNAs associated with biomineralization in gastropods, by far the largest class of molluscs.
Here we have surveyed the miRNA complement from a number of adult tissues in the freshwater pulmonate gastropod L. stagnalis. In addition we focus on miRNAs involved in biomineralization by inducing shell regeneration. Specifically, we compare the miRNA complements of shell-regenerating mantle regions against immediately adjacent, relatively inactive biomineralizing mantle regions. This matched pairs design allows us to identify miRNAs that potentially regulate either inhibitors or enhancers of shell formation. In this way we identify a number of novel and conserved miRNAs that are likely to be regulating shell forming genes. To support their validity as genuine miRNA molecules we characterize some of the primary genomic features and we assess their degree of evolutionary conservation by comparison to miRBase and mirGeneDB2.0 entries. To our knowledge this is the first report to identity miRNAs associated with the process of biomineralization in a gastropod.

3
under the International Union for Conservation of Nature (IUCN's) list of threatened species. As a non-cephalopod mollusc this work was exempt from regulations outlined by the University of Göttingen Ethics Committee. We also applied the 3R principles (replace, reduce, refine) in all of our animal work.
Adult Lymnaea stagnalis were maintained in our laboratory aquaria as previously described [40]. In order to initiate shell repair a notch was carefully ground into the growing edge of the right side of the shell of 10 individuals using a Dremel 8200. The notch and the immediately adjacent undamaged shell edge were then marked with a permanent marker in order to accurately monitor shell regeneration relative to "non-regenerating" adjacent shell. Once approximately 50% of the notch was deemed to have been repaired (for most snails this took 12-36 h), snails were euthanized by swift decapitation and bisection of the cerebral ganglia, and the mantle tissue immediately underlying the regenerating shell and underlying the adjacent non-regenerating shell were collected for RNA extraction. Note that we were mindful to only use snails that were efficiently repairing the notch and had not deposited any appreciable shell material in the adjacent shell edge (Fig. 1). In this way four biological replicates (i.e. four animals providing four matched pairs of regenerating and non-regenerating mantle samples) were obtained. In addition we performed small RNASeq on seven adult tissues (not replicated): mantle edge (distal); mantle proximal; cephalic tentacle; cephalic lobe; cephalic ganglia; foot; buccal mass. Total RNA was extracted from all samples using Qiazol (Qiagen) following the manufacturers protocol. TruSeq small RNA libraries were constructed and 50 bp single end sequencing was performed on the Illumina HiSeq2000 platform by the Kiel sequencing centre for the seven adult tissues and by the TAL (Göttingen) for the matched-pair mantle samples. Raw data are available from the NCBI SRA under the following accession numbers: SAMN15898436-SAMN15898450.

Sequence data processing
All Illumina raw data files were initially processed by the command line version of FastQC v.0.11.9 [41]. Remaining adaptor sequences were trimmed using BBDuk from the BBTools package [42] using the default Illumina adaptor file, a kmer of 25 and a minimum kmer of 8. Trimmomatic v0.36 was then used to quality trim all reads with leading and trailing values of 5, window size of 4 bases, quality of 15, target length of 40 bp, strictness of 0.75 and a minimum length of 15 bp. Reads were then error corrected with Tadpole from the BBTools package [42] using a kmer of 50 bp. Potentially contaminating rRNA reads were removed using BBDuk from the BBTools package (42) using a ribo-kmer library derived from the Silva database [43] and publicly provided by Brian Bushnell [44] with a kmer of 31 bp. Finally reads were filtered on length (minimum 18 bp, maximum 33 bp). All filtered reads were then processed again by FastQC v.0.11.9 [41] in order to assess the overall impact of these filtering measures. The detailed results of these filtering steps are available in Additional file 1: Table S1.

Identification of miRNA synthesis pathway genes
We used transcriptome data from one of our previous studies [40] to search for miRNA synthesis genes. We used 30 miRNA synthesis genes previously reported from the bivalve Mytilus galloprovincialis [18] as queries. These included Dicer (KT447258.1), Drosha (KT447251.1) and DGCR8, also known as Pasha (KT447252.1). All 30 protein sequences were compared against the L. stagnalis transcriptome assembly using tBLASTn implemented in the BLAST v2.2.28+ package.

Differential miRNA expression during shell repair, target prediction and ISH
Short RNA reads from the shell regeneration experiment were used for differential expression analysis using the CLC Genomics Workbench (version 12). Short reads were first mapped to the predicted miRNAs (see above) using the 'RNA-seq analysis' function. The mismatch cost was set to 2, insertion and deletion cost to 3, length and similarity fraction to 0.8 and the maximum number of hits for a single read was set to 30. The expression level of each miRNA was estimated from total read counts. Differential expression of miRNAs was performed using the "Small RNA" module and the trimmed mean of M (TMM) values was used for normalization [49]. All comparisons were paired (i.e. 1 animal provided 1 matched pair of samples) with the shell-regenerating mantle as the treatment and the immediately adjacent (non-shell-regenerating) mantle as the control. False discovery rate (FDR) corrected p-values were used for all statistical comparisons. In order to identify putative targets of all differentially expressed miRNAs we employed miRanda v3.3a [50]. We used a score threshold higher than 140 (-sc 80) and an energy threshold lower than − 20 kcal/mol (-en − 20). The "Gap Open Penalty" was fixed to −9, "Gap Extend Penalty" to −4 and the "scaling parameter" to 4. In situ hybridisation against tissue sections was performed as previously described [40].

Characterization of miRNA production machinery
Our preliminary survey of the miRNA synthesis machinery present in a master transcriptome dataset of L. stagnalis returned sequences with significant similarity to all 30 query sequences (Additional file 2: Table 2). As expected, this suggests that all of the components necessary for the generation of miRNAs are present in the L. stagnalis genome and are actively expressed. These sequences are available from GenBank under the following accession numbers: MT947760-MT947789.

miRNA RNA detection
A total of 102,578,980 short RNA reads derived from 7 adult tissues and 8 shell-regenerating samples (4 biological replicates of 2 matched paired conditions) passed our quality filters. These reads were used to identify potential miRNA precursors using miRDeep2 [45]. This process includes identifying the mature miRNA, Dicer cutting sites and star and loop sequences. With this stringent approach we identified 394 putative precursor miRNAs (pre-miRNAs) (Additional files 3, 4: Files 1 and 2). Of these 394 pre-miRNAs 370 are unique (i.e. 48 have identical sequences and have apparently been duplicated within the genome). These 370 unique pre-miRNAs generate 336 unique mature miRNAs (i.e. 68 unique premi-RNAs will generate a mature miRNA that has an identical sequence to another mature miRNA derived from another pre-miRNA; Additional file 5: File 3). The number of putative miRNAs we detect (394) can be compared to results recently reported for L. stagnalis (482) by Walker et al. [38], however only 81 of our miRNAs have a 100% to those reported by Walker et al. (Additional file 6: Table S3). We suspect that the more stringent criteria we applied to our raw data (Walker et al. did not use mirDeep2 to verify short read candidates against a genome) have revealed a more robust collection of miRNAs for L. stagnalis. Unfortunately we could not re-analyse the 482 miRNAs proposed by Walker et al. as their data is not publicly available. Nonetheless both of these miRNA complements for L. stagnalis are higher than the 140 reported for Caenhorbditis elegans [51], but lower than the 2,300 reported for humans [52]. We mapped the 102,578,980 filtered short RNA reads to these 394 reference miRNAs and found that 74.84% mapped. Many reads mapped to lst-mir-14694 (25,925,554 reads, accounting for 25% of all reads), indicating that relatively few miRNAs are responsible for the majority of miRNA expression. Of the short reads that could be mapped to their respective pre-miRNAs, the majority (91.13%) located to the mature region, while very few (0.61%) mapped to the loop, and the star region (8.27%) giving us further confidence that the majority of our raw data is biologically genuine. Among our 394 mature miRNAs, 110 (30%) share similarity with a miRBase entry: 7 have similarity with H. rufescens miRNAs, 82 with L. gigantea miRNAs and 21 with M. leonina miRNAs, leaving 284 (72.08%) novel and potentially unique to L. stagnalis. Of these 284 novel miRNAs 4 were previously identified by Walker et al. [38]. A similar search against the mirGeneDB2.0 database revealed that 96 of the 394 (24%) mature miRNAs have been identified in a diversity of other animals with most of the top hits against molluscs such as Lottia and Crassostrea, but also against vertebrates (Xenopus tropicalis) insects (Tribolium castaneum) and echinoderms (Strongylocentrotus purpuratus) (Additional file 6: Table S3). To put these numbers in context we also performed some additional searches of mirGeneDB2.0 using datasets from two molluscs and one brachiopod. For 81 miRNAs from the gastropod Lottia we could identify 69 hits (85%), for 155 miRNAs from the oyster Crassostrea we could identify 77 (50%) hits, and for 109 miRNAs from the brachiopod Lingula we could identify 97 (89%) hits. This suggests that Lymnaea has a relatively high proportion of evolutionarily novel miRNAs. Interestingly, the 70 most highly expressed mature miRNAs (those with more than 10,000 mapped reads) accounted for 99.77 % (76,593,842 reads) of all mapped reads. All 70 of these highly expressed miRNAs, except one, were present in all surveyed tissues, and the majority share sequence similarity with a miRBase v22 (66/70) or mirGeneDB2.0 (65/70) entry (Additional file 6: Table S3).

Differential expression of miRNAs involved in shell regeneration
In order to specifically focus on miRNAs regulating biomineralization, we searched for miRNAs differentially expressed in mantle tissue actively regenerating shell material. We identified 10 miRNAs that were significantly differentially expressed (FDR p-value < 0.05): 4 were down-regulated in shell-regenerating mantle tissue, while 6 were up-regulated in shellregenerating mantle tissue ( Table 1). Two of these have total read counts higher than 200,000 (Table 1). For 8 of these 10 differentially expressed miRNAs the majority (> 90 %) of reads that map to them are located in the mature region of the miRNA ( Table 1). Five of these 10 differentially expressed miRNAs also have similarity to a miRBase entry, and 9 to a mirGeneDB2.0 entry, suggesting they are evolutionarily conserved. Interestingly, none of these differentially expressed miRNAs have similarity with functionally characterized miRNAs, or with 12 miRNAs previously identified as playing a role in biomineralization in bivalves [24][25][26][27][28][29][30][31]. We also took these 12 miRNAs from previous reports focused exclusively on bivalves, and searched for potential orthologs in all of our miRNA data. Among these, 2 have perfect matches to miRNAs in our dataset (miR-1990c-3p and miR-876, [24]), one has a perfect match but is one base shorter (miR-29a [27]), and two have a single mismatch (miR-183, [30]; miR-9b-5p [28]). Finally, one of these 12 bivalve miRNAs (pm-miR-2386 [29]) is present in the L. stagnalis genome with 100 % identity, but was not present in any of our short read data.
We also surveyed our data for miRNAs that might play a conserved role in wound repair (although we were extremely careful not to damage the mantle tissue when applying the notch to the shell). Indeed we found none of the miRNAs differentially expressed during shell repair were associated with previously described roles in wound repair in other metazoans [53][54][55][56][57][58].

Genomic architecture of differentially expressed miRNAs
Among the 10 differentially expressed miRNAs, two are located within ~ 4,400 bp of each other. Lst-miR-4873 and lst-miR-4881 are located on either side but on the non-coding strand of a short (246 bp) single-exon predicted coding sequence (CDS; Fig. 2). Lst-miR-4873 is localised 852 bp upstream and lst-miR-4881 3322 bp downstream of the CDS. Both miRNAs have similar read counts (2725 reads for lst-miR-4873 and 2366 for lst-miR-4881 ( may be expressed in a single transcriptional unit. The CDS encodes an apparently novel protein (it contains no recognizable domains, and does not share any similarity with sequences in the non-redundant NCBI database). Four of the 10 differentially expressed miRNAs are located within exons: lst-miR-3861, lst-miR-16924, lst-miR 26499 and lst-miR-33878 (Fig. 2). Of these 4 miRNAs only lst-miR-3861 is located on the non-coding strand. Lst-miR-18140 is located within an intron on the coding strand of a predicted novel gene, while Lst-miR-8292 is also located within an intron of a predicted novel gene but on the non-coding strand (Fig. 2). Lst-miR-3847 is located 895 bp upstream on the non-coding strand of another novel gene, and finally lst-miR-28733 is located in a genomic area where no CDS was detected (Fig. 2).

Predicted targets of ten differentially expressed miRNAs
In order to identify potential targets of the ten differentially regulated miRNAs involved in shell repair and biomineralization we employed miRanda [50] and supplied it with a set of 181,393 annotated genes from the L. stagnalis genome. This gene set included our 32 previously in situ verified L. stagnalis shell forming genes [40] and a set of 32 genes associated with left/right shell asymmetry in L. stagnalis [59]. While none of the top putative targets could clearly be assigned a functional role in shell repair or biomineralization (Table 1), the most down-regulated miRNA (lst-miR-4873 at −3.4 fold) was predicted to target Dermatopontin as the fourth top target (Additional file 7: Table S4). Dermatopontin was first reported as a major protein from the shell of the freshwater gastropod Biomphalaria glabrata [60] and has since been identified in other gastropods [61] and bivalves [62]. Given that lst-mir-4873 is significantly downregulated in mantle tissue that is actively repairing shell, it would follow that targets of lst-mir-4873 would therefore be significantly more abundant. This prompted us to perform in situ hybridization on mantle tissue of healthy juveniles against Lst-dermatopontin. As expected, we found that indeed this gene has a specific and strong spatial expression pattern in the belt region of the mantle (Fig. 3). While this was the only clear relationship between a significantly differentially regulated miRNA and a functionally well-defined shell matrix protein (SMP) many of the other targets listed in Additional file 7: Table S4 have no sequence similarity with Swissprot entries, a trait common to SMPs. Assessing the spatial expression profiles of all 200 of these putative targets would be a considerable undertaking, however this approach may reveal additional components of the biomineralizing toolkit.

Conclusions
We have identified what is likely to be a significant proportion of the actively expressed miRNAs in adult tissues form the Great Pond Snail Lymnaea stagnalis. A large fraction of these (more than 70 %) are apparently novel with no clear homologs in miRBase or MirGeneDB2.0. Through the application of a shell-regeneration assay we have identified 10 miRNAs that are associated with regulating the process of shell repair and biomineralization. The majority of the putative Fig. 3 In situ hybridization spatially localizes a shell-forming target of lst-mir-4873, Lstdermatopontin, to the mantle. The dashed-boxed indicates the magnified area shown in the inset and consists of the characteristic 'belt' region where many shell forming genes are expressed in the tall columnar cells. The inset image has also been stained with DAPI to reveal cell nuclei (light blue), while the asterisk indicates the periostracal groove