1. Introduction {#s0005}
===============

MicroRNAs (miRNAs) are a class of about 18--24 nucleotide non-coding RNA molecules that are negative regulators of target genes via complementary base-pair interactions in a variety of eukaryotic organisms ([@b0140]). MiRNAs originate from the primary transcripts (pri-miRNA), which were transcribed by polymerase II or III and then processed by RNAase-III-like enzyme (Drosha) to produce hairpin precursor miRNAs (pre-miRNA). Pre-miRNAs are then exported into the cytoplasm, wherein pre-miRNAs are further cleaved by a second RNAase-III-like enzyme (Dicer) to form mature miRNA ([@b0050]).

In order to identify newer miRNAs, direct cloning, solexa deep sequencing and computational techniques have been widely used in many species. The direct cloning method is considered as an efficient way to obtain highly expressed miRNAs, while it is difficult to find miRNAs expressed at a low level. Solexa deep sequencing is preferably used to find lowly expressed miRNAs and species-specific miRNAs, while some miRNAs are hard to sequence due to their physical properties or post-transcriptional modifications, such as editing or methylation. It is well known that many miRNAs are evolutionarily conserved from species to species. This suggests a powerful strategy for identifying miRNAs in specific species through the computational technique based on the highly conserved sequences in mature miRNAs. And many miRNAs in a variety of species have been successfully identified by computational searches ([@b0040], [@b0055], [@b0115], [@b0150]).

*Pinctada martensii*, one of the most important components of molluscan mariculture, was primarily cultured for pearl production in China and Japan. It is an outstanding model organism for elucidating the mechanism underlying biomineralization. In our previous research, using deep solexa sequencing technology, 258 *P. martensii* miRNAs (pm-miRNA) were identified ([@b0065]). However, the number of the identified pm-miRNA was still less than in other animal species. In this study, we hope to get more candidate pm-miRNAs using computational technology based on the transcriptome database of pearl sac from *P. martensii* in our previous research ([@b0160]). Furthermore, we assessed the potential function of the predicted pm-miRNAs in biomineralization. This study provided insights into the mechanism underlying biomineralization in pearl and shell formation of pearl oyster *P*. *martensii*.

2. Materials and methods {#s0010}
========================

2.1. Experimental samples and RNA extraction {#s0015}
--------------------------------------------

Pearl oyster *P. martensii* (about 2 years of age) were collected from Liushagang, Zhanjiang, Guangdong Province of China. The pearl sacs were collected from pearl oyster and immediately stored in liquid nitrogen until used. Total RNAs were extracted using Trizol reagent. The integrity of RNA was determined by 1.2% formaldehyde-denatured agarose gel and staining with ethylene bromide. The quantity of RNA was determined by measuring OD260/OD280 with NanoDrop ND1000 Spectrophotometer.

2.2. Reference sequence data and employed software {#s0020}
--------------------------------------------------

The animal mature miRNAs were downloaded from miRNA database (<http://www.mirbase.org/>, released 20). Biomineralization-related genes were downloaded from NCBI GenBank. Comparative software (BLASTN) was downloaded from NCBI GenBank. The secondary structures of miRNA precursors were processed by Mfold ([http://mfold.rna.albany.edu/?q=mfold/RNA-Folding-Form](http://www.mfold.rna.albany.edu/?q=mfold/RNA-Folding-Form){#ir0035}). Target analyses between miRNAs and biomineralization-related genes were analyzed by microTar ([http://tiger.dbs.nus.edu.sg/microtar /running.html](http://tiger.dbs.nus.edu.sg/microtar%20/running.html){#ir0015}), miRanda algorithm (<http://www.microrna.org/microrna/getDownloads.do>) and RNA22 (<https://cm.jefferson.edu/rna22v2/>).

2.3. Prediction of miRNAs {#s0025}
-------------------------

The sequences used for pm-miRNA prediction were derived from the transcriptome database of pearl sac from *P. martensii* ([@b0160]). Procedures of searching for candidate pm-miRNAs were performed as described previously with some modification. Five criteria were used to predict mature pm-miRNAs and pre-pm-miRNAs: (1) predicted mature pm-miRNAs were allowed to have only 0--4 bp mismatches in sequence with known animal mature miRNAs; (2) the mismatched nucleotides were not permitted in the miRNA seed sequence (2--8 bp); (3) pre-pm-miRNAs sequence can fold into an appropriate hairpin secondary structure that contains mature pm-miRNA sequence within one arm of the hairpin and have the smallest folding energy; (4) the predicted pm-miRNA hairpin is on the same stem of the hairpin as the blast source of the candidate miRNA hairpin; (5) miRNA::miRNA^∗^ duplex mismatches were restricted to eight or fewer.

2.4. Identification of pm-miRNAs using stem-loop qRT-PCR {#s0030}
--------------------------------------------------------

Stem loop qRT-PCR analysis was employed to validate and determine the specific expression of pm-miRNAs with 5.8S as the internal control ([@b0070]). Stem-loop RT primers, real-time PCR primers were designed as previously described ([@b0010]). All the primers used in this study are listed in [Table 1](#t0005){ref-type="table"}. The relative expression of mature pm-miRNAs was calculated with the 2^(Δct)^, and the relative expression of primary miRNAs was represented with reads per kilo bases per million reads (RPKM) calculated by reads in pearl sac transcriptome.

2.5. Target prediction for biomineralization-related genes {#s0035}
----------------------------------------------------------

Biomineralization-related genes in *P. martensii* were downloaded from NCBI GenBank. Target analysis between the candidate pm-miRNAs and the biomineralization-related genes was performed using microTar ([@b0130]), miRanda ([@b0025]) and RNA22 ([@b0090]). When used microtar algorithm, the dimer--monomer energy difference (*g*) was used to estimate the necessary in duplex formation. It was considered as the functional targets only when *g* ⩽ −5. Functional target genes of miRNAs using the miRanda algorithm were restricted with score ⩾100 and free energy ⩽−10 kcal/mol. And the RNA22 was used to determine the most favorable hybridization site between miRNA and mRNA with sensitivity of 63%, specificity of 61%, and seed size of 7, allowing the maximum of 1 un-paired bases in the seed, the minimum number of paired-up bases in heteroduplex ⩾12, the maximum folding energy for heteroduplex ⩽−12 kcal/mol.

3. Results {#s0040}
==========

3.1. Identification of candidate pm-miRNAs {#s0045}
------------------------------------------

After eliminating the reported pm-miRNAs, a total of 30 pm-miRNAs were identified according to the criteria described in Section [2](#s0010){ref-type="sec"}. The length range of the predicted pm-miRNAs was 17--25 bp. The number of sequences with 19 bp and 22 bp was significantly higher than others ([Fig. 1](#f0005){ref-type="fig"}). The free energies of the secondary hairpin structures ranged from −13.20 to −30.0 kcal/mol. The hairpin structures of those are shown in [Fig. S1](#s0075){ref-type="sec"}. Among the identified 30 pm-miRNAs, 13 pm-miRNAs were generated from the complementary strand of protein-coding mRNAs by blastx analysis ([Table S1](#s0075){ref-type="sec"}), indicating the potential origination of some pm-miRNAs. And the other 17 pm-miRNAs could not be annotated using blastx and tblastn analysis. Notably, 10 of the 30 pm-miRNAs were homologous with those reported in our previous research. To discriminate these pm-miRNAs from the reported ones, we designated these pm-miRNAs as pm-miR-1b, pm-miR-13b and pm-miR-138b, pm-miR-205b, pm-miR-1638b, pm-miR-989b, pm-miR-4039b, pm-miR-135b, pm-miR-375b and pm-miR-1892b. The remaining 20 miRNAs were novel pm-miRNAs without identity to the reported pm-miRNAs.

3.2. Validation of pm-miRNAs by stem-loop qRT-PCR {#s0050}
-------------------------------------------------

In order to validate the presence of the candidate pm-miRNAs in this study, we randomly selected eight pm-miRNAs and performed stem loop qRT-PCR analyses using 5.8S as the internal reference gene. All of the selected pm-miRNAs could be detected easily with the maximal CT = 30. Pm-miR-1b (CT = 13) and pm-miR-669n (CT = 14) showed the highest expression levels compared with other pm-miRNAs.

All the obtained pm-miRNAs were predicted from the unigenes in the transcriptome database of pearl sac in *P. martensii*, thus the RPKM of the unigenes represented the expression level of pri-miRNAs. To investigate the expression correlation between the pri-pm-miRNA and mature pm-miRNA, we calculated the proportion of each pri-pm-miRNA or mature pm-miRNA in the correspondent aggregate sum of pri-pm-miRNAs or mature pm-miRNAs, respectively. Results showed only two relatively highly expressed miRNA, pm-miR-135 and pm-miR-748, had some expression correlation between the pri-pm-miRNA and mature pm-miRNA ([Fig. 2](#f0010){ref-type="fig"}). This result suggested that the consistency of the expression profile between pri-miRNA and mature miRNA might be presented in some highly expressed miRNAs, an indicator of a complex process of miRNA accumulation and degradation.

3.3. Functional prediction of pm-miRNAs in biomineralization {#s0055}
------------------------------------------------------------

To elucidate the function of the obtained pm-miRNAs in biomineralization, the putative target associations between pm-miRNAs and the reported biomineralization related-genes were analyzed by microtar, mirRanda and RNA22. Sixty biomineralization related-genes with completed 3′ UTR sequences were downloaded from NCBI database. All of three programs were used to predict the miRNA-mRNA interaction sites. As shown in [Table S2](#s0075){ref-type="sec"}, hundreds of miRNA--mRNA interaction sites were obtained by only one program, respectively. And the interaction sites vary widely based on different prediction programs. For example, by miRanda, N14 family genes were predicted to be regulated by pm-miR-2386 with the same target site (GAACGAAAC, [Fig. 3](#f0015){ref-type="fig"}a), while, by microTar, five of the six N14 family genes were potentially regulated by pm-miR-1892 with the common target site (UCUCAGA, [Fig. 3](#f0015){ref-type="fig"}b) and low dimer-monomer energy difference (no more than −5). Both results seemed interesting while they need further experimental elucidation.

Among these interaction sites, only thirteen biomineralization-related genes that are potentially regulated by multiple pm-miRNAs were predicted by three tools simultaneously ([Table S2](#s0075){ref-type="sec"}). Even so, a complex multi-gene regulation system existed between pm-miRNAs and the biomineralization related-genes. For example, pm-miR-2386 was indicated to regulate five genes, such as pearlin, N16-7 and chitin synthase. The Glycine-rich protein named insoluble protein (ISP), one molluscan shell framework proteins, could be regulated by two miRNAs, pm-miR-1719 and pm-miR-669n ([Fig. 4](#f0020){ref-type="fig"}).

In addition, analyzing these interactions, we found the interactive sites predicted by different programs even between the same miRNA--mRNA, such as pm-miR-13a-BMP2, were different, as shown in [Fig. 5](#f0025){ref-type="fig"}, indicating the striking difference in algorithm of three programs.

4. Discussion {#s0060}
=============

Increasing evidence has demonstrated that miRNAs participate in various biological processes by negatively regulating the target mRNA. As their prominent functions, hundreds of miRNAs have been identified in recent years, while only a very small amount of miRNAs in invertebrate especially in pearl oyster has been discovered and functionally identified. Computational strategy for discovery of miRNA is an efficient method that identifies miRNA in species without genome sequence and has been successfully applied to many species, such as insects, elegance and humans ([@b0020], [@b0035], [@b0045], [@b0075]). Using deep solexa sequencing, 258 pm-miRNAs were identified in our previous research ([@b0065]). In this study, a total of 30 candidate pm-miRNAs were identified using computational technology based on the transcriptome database of pearl sac from *P. martensii* ([@b0160]).

To validate the existence of the predicted pm-miRNAs, we performed stem-loop qRT-PCR, which is a reliable method to detect and measure the expression levels of pm-miRNAs. All the selected pm-miRNAs could be easily detected. Pm-miR-1b was presented with the most expression level, similar to the tested pm-miR-1a by deep sequencing in our previous research ([@b0065]). In mammals, miR-1 was specifically-expressed in muscle and involved in regulating muscle proliferation and differentiation ([@b0015]). It was reported that miR-1 could regulate calcium signaling pathways by targeting calmodulin ([@b0060]). Here, we also found interaction between pm-miR-1b and calmodulin by microTar program, indicating its functions in calcium signaling pathway in *P. martensii*.

We also compared the expression correlation between the pri-pm-miRNA and mature pm-miRNA. Results showed only two relatively highly expressed pm-miRNAs showed some consistency between primary and mature type. As is illustrated in the introduction, the maturation of miRNA required to undergo the process from the pri-miRNA to precursor by Drosha, and then by Dicer to the mature miRNA. And stability and decay of miRNAs were also affected by cell-cycle regulation, target binding and uridylation modifications ([@b0155]). These intricate processes may be responsible for the inconsistency between primary and mature miRNAs. Therefore, the expression level from RNA-seq representing primary miRNAs cannot display the final expression of mature miRNAs. The regulation of miRNA accumulation and decay needs to be unveiled.

In biominerals, the small amount of organic components, including polysaccharides, proteins, glycoproteins and proteoglycans, exerts crucial rules in crystal morphology, crystallographic orientation, and the superior material property of the formed crystal ([@b0085]). Meanwhile, biomineralization is a very complex and precise process in which the expression of each related protein is subject to fine regulation. Consequently, the regulators involved in biomineralization should also be paid attention. Generally, the function of miRNA is mediated by silencing gene expression of the targeted mRNAs. Target prediction tools were used to predict the functions of miRNAs and had been widely used in many reports especially in the initial research when you have little information about the miRNAs in this species ([@b0105]). To illustrate the potential function of the pm-miRNAs in biomineralization, we analyzed the interaction site between the obtained pm-miRNAs and the reported biomineralization-related genes by target prediction programs. As the difference in algorithm of different programs, the predicted interaction sites by different programs vary widely, which could be seen in our results. To reduce the false positive site, we screened the interaction site by three programs, including microTar, miRanda and RNA22. All of those tools were not relying on evolutionary conservation to discern functional targets, which made it possible to predict both conserved and non-conserved targets.

A total of thirteen biomineralization-related genes, predicted by three tools simultaneously, were potentially regulated by multiple pm-miRNAs. Both of pearlin and N16, belonging to the N16 family ([@b0030]), were matrix proteins isolated from the nacre of the molluscan shell. Pearlin, with a high proportion of Gly, Tyr, Asn and Cys, was presumed to bind to calcium as a template for nucleation during nacre formation ([@b0110]). In vitro studies indicated N16 could induce the formation of aragonite crystals ([@b0100]). Our results indicated pearlin and N16-7 were both potentially targeted by pm-miR-2386, implying a critical role of miR-2386 in the regulation of biomineralization. Previous studies showed that chitin was not only one of the important components of the periostracum but also one abundant organic matrix in the nacre and prismatic layer ([@b0080], [@b0120]). The chitin synthase, as the vital enzyme to synthesize the chitin ([@b0120], [@b0135]), was also inferred to be targeted by pm-miR-2386. In addition, the tyrosinase-like protein 1, as a member of the tyrosinase family which has been declared to be potentially essential for shell matrix maturation ([@b0005], [@b0145], [@b0150], [@b0155]), was also predicted to be negatively regulated by pm-miR-2386. Therefore, we presumed that pm-miR-2386 may be as the hub gene in the regulation network during shell formation by directly and indirectly controlling the organic matrix formation.

Generally speaking, matrix proteins were secreted from the epidermic cells of mantle tissue in the shelled mollusk. And similar to vertebrates, cell proliferation and differentiation in invertebrates also need the regulation of growth factors. Previous studies showed that bone morphogenetic protein 2 (BMP2), belonging to the transforming growth factor type beta supergene family, may have a key role in nacre formation. It was pointed out that the function of BMP2 in the formation of hard tissues was conserved from species to species ([@b0095], [@b0125]). Our results showed BMP2 was predicted to be regulated by pm-miR-13b. These results suggested that pm-miRNAs were involved in the biomineralization process not only through targeting the biomineralization-related proteins, but also through the controlling of cell proliferation and differentiation of biomineralization-related cells.

5. Conclusion {#s0065}
=============

A total of 30 pm-miRNAs and their precursors were identified based on the transcriptome data of *P. martensii*. Eight randomly selected pm-miRNAs were confirmed by stem-loop qRT-PCR. Some of the pm-miRNAs, such as pm-miR-2386 and pm-miR-13b, were predicted to participate in biomineralization by regulating the expression of biomineralization-related genes and the cell proliferation and differentiation of biomineralization-related cells. Therefore, this study enriched the miRNA databases of pearl oyster and provided an overview of the conserved pm-miRNAs as well as their potential functions in biomineralization.

Appendix A. Supplementary data {#s0075}
==============================

Supplementary Fig. S1Hairpin structures of the candidate pm-miRNAs predicted by bioinformatics analysis. The regions for generating mature miRNAs were marked by lines.Supplementary Table S1Candidate miRNAs predicted by computational identification.Supplementary Table S2pm-miRNAs functioning in biomineralization predicted by microTar miRanda and RNA22.

The studies were financially supported by grants of the National Natural Science Foundation of China (41206141, 31272635 and 31372526), Natural Foundation of Guangdong Ocean University (1212318), Excellent Young Scientist Foundation of Guangdong Ocean University (201402).

Peer review under responsibility of King Saud University.

Supplementary data associated with this article can be found, in the online version, at [http://dx.doi.org/10.1016/j.sjbs.2015.04.001](10.1016/j.sjbs.2015.04.001){#ir0030}.

![Size distribution pattern of the predicted pm-miRNAs.](gr1){#f0005}

![Expression Levels of eight randomly selected pm-miRNAs. Mature pm-miRNAs were detected by stem loop qRT-PCR using 5.8S as the internal reference gene. Pri-pm-miRNAs were validated by the RPKM of the unigenes in the transcriptome database of pearl sac in *P. martensii*.](gr2){#f0010}

![The target sites of pm-miR-2386 and pm-miR-1892b in the N14 family genes. The target sites of pm-miR-2386 (a) and pm-miR-1892b (b), indicated by black box, in the 3′ UTR regions of six of N14 family genes were predicted by miRanda and microTar, respectively.](gr3){#f0015}

![The target genes of pm-miR-2386, pm-miR-1719 and mir-669n. These target interactions were predicted by three programs, including microTar, miRanda and RNA22, simultaneously. TYR-like protein 1 represents tyrosinase-like protein 1; PfCHS1 represents the chitin synthase 1 of *P. fucata*; MSI60RP represents MSI60-related protein; ISP represents the gene named insoluble protein.](gr4){#f0020}

![Hairpin structure of pm-miR-13b and the potential target sites between pm-miR-13b and BMP2. (a) Hairpin structure of pm-miR-13b was predicted by M-fold. Target sites were predicted by microTar (b), miRanda (c) and RNA22 (d) programs.](gr5){#f0025}

###### 

Sequences of primers for stem-loop RT-PCR.

  Name          RT primer                                                    qPCR-antisense primer   qPCR-sense primer
  ------------- ------------------------------------------------------------ ----------------------- ------------------------
  pm-miR-765    GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACCACCTCCT     TGCGTGTCGTGGAGTC        GGAGGAGAAGGAGGAG
  pm-miR-669n   GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACACACACAT     TGCGTGTCGTGGAGTC        ATTTGTGTGTGTTATGTGTGT
  pm-miR-205b   GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACCAGAGGCC     TGCGTGTCGTGGAGTC        ACCTTCATTCCACCGGCCTC
  pm-miR-135b   GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACCACACAGG     TGCGTGTCGTGGAGTC        TATGGCTTTTCTCCTGTGTG
  pm-miR-748    GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACCTCTTCAC     TGCGTGTCGTGGAGTC        TGGACGGAAGTGTGAAGAG
  pm-miR-3377   GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACTATACAGA     TGCGTGTCGTGGAGTC        TGAAGTTATGACGTCTGTATA
  pm-miR-3122   GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACAGACTGTC     TGCGTGTCGTGGAGTC        GTTGGGACAAGAATGGACAGT
  pm-mir-1b     GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACTCCATACT     TGCGTGTCGTGGAGTC        TGGGATAGCAAAGAAGTATGGA
  5.8S          GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACACCGACCCTG   TGCGTGTCGTGGAGTC        GCTAGCTGCGTGAATTAATGTG

[^1]: These authors contribute equally.
