Computational prediction of candidate miRNAs and their potential functions in biomineralization in pearl oyster Pinctada martensii.

MicroRNAs (miRNAs) are a class of non-coding RNA molecules with presumed post-transcriptional regulatory activity in various biological processes, such as development and biomineralization. Pinctada martensii is one of the main species cultured for marine pearl production in China and Japan. In our previous research, 258 pm-miRNAs had been identified by solexa deep sequencing in P. martensii, while it is far from the number of miRNAs found in other species. In this study, based on the transcriptome database of pearl sac, we identified 30 candidate pm-miRNAs by computational prediction. Among the obtained 30 pm-miRNAs, 13 pm-miRNAs were generated from the complementary strand of protein-coding mRNAs, and 17 pm-miRNAs could not be annotated using blastx and tblastn analysis. Notably, 10 of the 30 pm-miRNAs, such as pm-miR-1b, pm-miR-205b and pm-miR-375b, were homologous with the reported pm-miRNAs, respectively. To validate the existence of the identified pm-miRNAs, eight randomly selected pm-miRNAs were tested by stem loop quantitative RT-PCR analyses using 5.8S as the internal reference gene. Target prediction between the obtained pm-miRNAs and biomineralization-related genes by microTar, miRanda and RNA22 indicated pm-miR-2386 and pm-miR-13b may be the key factors in the regulation network by regulating the formation of organic matrix or the differentiation of mineralogenic cell during shell formation. Thus, this study enriched miRNA databases of pearl oyster and provided a new way to understand biomineralization.

Introduction

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 (Yu et al., 2012). 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 (Hinske et al., 2010).

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 (He et al., 2008a, Huang et al., 2010, Sheng et al., 2011, Zhang et al., 2012b).

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 (Jiao et al., 2014). 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 (Zhao et al., 2012). 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.

Materials and methods

Experimental samples and RNA extraction

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.

Reference sequence data and employed software

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). Target analyses between miRNAs and biomineralization-related genes were analyzed by microTar (http://tiger.dbs.nus.edu.sg/microtar /running.html), miRanda algorithm (http://www.microrna.org/microrna/getDownloads.do) and RNA22 (https://cm.jefferson.edu/rna22v2/).

Prediction of miRNAs

The sequences used for pm-miRNA prediction were derived from the transcriptome database of pearl sac from P. martensii (Zhao et al., 2012). 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.

Identification of pm-miRNAs using stem-loop qRT-PCR

Stem loop qRT-PCR analysis was employed to validate and determine the specific expression of pm-miRNAs with 5.8S as the internal control (Kolachala et al., 2010). Stem-loop RT primers, real-time PCR primers were designed as previously described (Chen et al., 2005a). All the primers used in this study are listed in Table 1. 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.

Table 1

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

Target prediction for biomineralization-related genes

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 (Thadani and Tammi, 2006), miRanda (Enright et al., 2004) and RNA22 (Miranda et al., 2006). 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.

Results

Identification of candidate pm-miRNAs

After eliminating the reported pm-miRNAs, a total of 30 pm-miRNAs were identified according to the criteria described in Section 2. 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). 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. Among the identified 30 pm-miRNAs, 13 pm-miRNAs were generated from the complementary strand of protein-coding mRNAs by blastx analysis (Table S1), 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.

Figure 1

Size distribution pattern of the predicted pm-miRNAs.

Validation of pm-miRNAs by stem-loop qRT-PCR

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). 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.

Figure 2

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.

Functional prediction of pm-miRNAs in biomineralization

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, 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. 3a), while, by microTar, five of the six N14 family genes were potentially regulated by pm-miR-1892 with the common target site (UCUCAGA, Fig. 3b) and low dimer-monomer energy difference (no more than −5). Both results seemed interesting while they need further experimental elucidation.

Figure 3

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.

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). 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).

Figure 4

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.

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, indicating the striking difference in algorithm of three programs.

Figure 5

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.

Discussion

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 (De Souza Gomes et al., 2013, Grad et al., 2003, He et al., 2008b, Lai et al., 2003). Using deep solexa sequencing, 258 pm-miRNAs were identified in our previous research (Jiao et al., 2014). 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 (Zhao et al., 2012).

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 (Jiao et al., 2014). In mammals, miR-1 was specifically-expressed in muscle and involved in regulating muscle proliferation and differentiation (Chen et al., 2005b). It was reported that miR-1 could regulate calcium signaling pathways by targeting calmodulin (Ikeda et al., 2009). 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 (Zhang et al., 2012c). 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 (Marin et al., 2012). 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 (Rajewsky, 2006). 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 (Fang et al., 2011), 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 (Samata et al., 1999). In vitro studies indicated N16 could induce the formation of aragonite crystals (Miyashita et al., 2000). 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 (Levi-Kalisman et al., 2001, Suzuki et al., 2007). The chitin synthase, as the vital enzyme to synthesize the chitin (Suzuki et al., 2007, Weiss et al., 2013), 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 (Aguilera et al., 2014, Zhang et al., 2012a, Zhang et al., 2012b, Zhang et al., 2012c), 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 (Miyashita et al., 2008, Takami et al., 2013). 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.

Conclusion

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.