1 INTRODUCTION

MicroRNAs are ~22 nucleotides (nt), single-stranded, endogenous and non-coding RNAs that play important roles in many physiological and developmental processes (Ambros, 2004; Bartel, 2004; He and Hannon, 2004). Since the first miRNA was discovered in the early 1990s (Lee et al., 1993), increasing evidence on multicellular eukaryotes has shown that miRNAs regulate basic cellular functions, such as proliferation, differentiation, and death (Hannon, 2002; Hwang and Mendell, 2006; Shivdasani, 2006; Lenkala et al., 2014; Yao, 2016), and have multipurpose biological functions, including development, growth, signalling, differentiation, protein degradation (Behm-Ansmant et al., 2006), stress response, disease development, and acting as a signature for cancer (Fleming et al., 2013; Tüfekci et al., 2014). miRNAs can negatively mediate gene expression by recognizing completely or partially complementary sequences in target genes. Previous studies show that miRNA genes are present in the introns of coding genes or intergenic regions. A series of steps was involved in the biogenesis of miRNAs. In the nucleus, miRNA genes are transcribed into several hundred nucleotide-long primary miRNAs (pri-miRNAs). Pri-miRNAs can liberate pre-miRNA (precursor miRNA) hairpin structures, which are then exported from the nucleus to the cytoplasm. Next, a pre-miRNA can produce a duplex intermediate (miRNA:miRNA^*) in the cytoplasm. Finally, mature miRNAs that usually come from one strand of the miRNA:miRNA^* duplex combine with Dicer enzyme and other associated proteins to form an active RNA-induced silencing complex (RISC). The RISC can repress translation by binding to the 3′UTR of the target mRNA (Hannon, 2002; Bartel, 2004; Yang and Wang, 2011).

As reported in previous studies, there are two prevalent methods for miRNA identification, the computational and sequencing approaches. Although miRNA sequencing provides more complete miRNA information, the computational method remains a prevalent strategy for miRNA prediction before high-throughput sequencing experiments are performed, especially for organisms with abundant genome and transcriptome resources. (Tong et al., 2006; Baev et al., 2009; Zhou et al., 2009; Frazier et al., 2010; Huang et al., 2010; Patanun et al., 2013; Akter et al., 2014; Han et al., 2014a, b; Li et al., 2016a). The homology-based and ab initio approaches are two main computational methods for miRNA identification. The homology-based approach takes advantage of the high conservation of miRNA among species. Although homology-based methods are convenient, they can only identify known miRNA molecules that have been previously reported. The ab initio approach relies on the sequence, genomic organization and structure features of miRNA precursors, which can detect miRNAs that lack obvious homology with those in other organisms.

Large yellow croaker (Larimichthys crocea), belonging to the order Perciformes, is a widely cultured and important economic marine fish in China (Zheng et al., 2006). Although miRNAs molecules were predicted for large yellow croaker (Huang et al., 2016), a previous study only identified miRNAs from the genome sequence. However, scientists have suggested that the transcriptome is a more reliable resource to identify expressed non-coding RNA (Wan et al., 2012; Ou et al., 2013; Xu et al., 2013, 2015; Prakash et al., 2016). Therefore, integrated miRNA identification by combining homology-based and ab initio methods from both the genome and transcriptome are needed for systematic analysis of miRNA in large yellow croaker.

A single nucleotide polymorphism (SNP) is an important small genetic variant in the population and provides valuable genetic resources for bio-diversity, genome evolution, and trait association analyses (Shastry, 2009; Nicoloso et al., 2010; Gong et al., 2012; Jiang et al., 2016; Li et al., 2016b). Many studies have revealed that a genomic SNP can affect many essential biological processes via gene expression regulation by miRNA-related interaction networks, especially for SNPs that reside in miRNAs and the binding region of mRNA. For example, rs2910164 of miR-146a plays an important role in gastric cancer (Jiang et al., 2016), and genetic variants in mir-499 and mir-608 could significantly increase the risk of lung cancer (Li et al., 2016b). Interestingly, miRNA seed regions, 2–8 nucleotides of the 5′ end, bind to complementary sites in the 3′-untranslated regions (3′UTRs) of target messenger RNAs (mRNAs) and play an important role in regulating expression during translation. Therefore, alterations in the miRNA seed sequences might dramatically influence the regulation network of miRNA. Previous studies have revealed that polymorphisms are rare in mature sequences, especially in seed regions because of the conservation of the mature miRNA sequence and its seed regions (Lewis et al., 2005). However, most studies were performed on human species, and little is known about the distribution of SNPs and their effects around miRNA in large yellow croaker and other teleost species.

In this study, an integrated computational method combining homology-based and ab initio predictions was used to identify potential miRNAs in large yellow croaker from both genome and transcriptome sequences. Therefore, conserved miRNAs and large yellow croaker-specific miRNAs were identified through this comprehensive method. For a better understanding of the miRNA regulation function, target genes of those miRNAs were predicted and annotated with the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. The distribution of SNP markers was investigated around miRNA sequences to understand how genetic variants influence the biological function of miRNAs in living cells. The SNP loci residing in miRNA and their mRNA binding sites were analysed to reveal their possible biological effects. This study provides a systematic method of miRNA detection and miRNA-related SNP analysis from the genome and transcriptome of large yellow croaker using the integrated strategy, providing a valuable reference for miRNA identification and functional investigation of other teleost species.

2 MATERIAL AND METHOD 2.1 Sequence data and reference miRNAs

The latest miRBase release (Release 21), containing 35 828 miRNA genes from 223 species, was downloaded from the miRBase (Kozomara and Griffiths-Jones, 2014), and a total of 26 332 animal miRNAs were retrieved as the reference set of the miRNA sequence. To omit redundant or overlapping miRNAs, repeated sequences of miRNAs were removed and the remaining sequences were used as query sequences for homology-based miRNA identification. The genome and transcriptome sequences of large yellow croaker were downloaded from the NCBI database with genebank accession number JRPU00000000 (Ao et al., 2015) and project accession number PRJNA254539 (Xiao et al., 2015), respectively. The transcriptome sequence came from various developmental stage, including embryos cells, larval, 11 juvenile and 2 adult (one male and one female), and different tissues including ocular, skin, muscle, gonadal, intestinal, liver, kidney, blood gall and air bladder.

2.2 The softwares used in the study

The alignment tool BLAST version 2.2.31 was used for conserved miRNA prediction. Secondary structures of pre-miRNAs were evaluated by RNAfold version 2.1.9 (http://www.tbi.univie.ac.at/RNA). RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de) and Miranda (http://www.microrna.org) were used for miRNA target gene prediction.

2.3 Identification of the conserved miRNAs

The workflow for identifying conserved miRNAs is presented in Fig. 1. The genome sequences of large yellow croaker were used for the homology search against the known unique mature animal miRNA sequences. The homology search was performed using BLASTN with an E-value of 10 and word-match size of 7. The following seven criteria were also considered to identify the candidate miRNAs from the alignment results:

1) the number of mismatches and gaps between known miRNAs and genomic sequences was less than 3;

2) the minimal length of the pre-miRNA was set as 55 nt;

3) the pre-miRNA was folded into an appropriate stem loop hairpin secondary structure;

4) the mature miRNA sequence and its opposite miRNA strand had less than 5 nt mis-matches;

5) the mature miRNA sequence resided in one arm of the hairpin structure;

6) the A+U content was between 30 and 70%;

7) both minimal folding free energy (MFE) and a minimal folding free energy index (MFEI) are essential to distinguish between miRNAs and other small RNAs. MFE and MFEI were lower than -20 and 0.8 kcal/mol, respectively.

To reduce the false positive rate for miRNA identification, the pre-miRNA sequences were blasted against the protein database with BLASTX to remove the coding sequence. The word match size between candidate pre-miRNA and protein sequences was kept at 7. For a more accurate result, the E-values of all of the alignment were collected to determine the best E-value using ROC analysis. A Support Vector Machine (SVM), exhibiting excellent pre-miRNA classification in previous studies for other species (Ng and Mishra 2007; Batuwita and Palade 2009; Ding et al., 2010), was used to scan the pre-miRNA sequences to distinguish true pre-miRNAs from false ones. Thirty-six pre-miRNA sequence features (see Supplementary Table 1 for detailed information) of 328 Danio rerio pre-miRNAs from miRBase and 350 random pre-miRNA sequences were used to train the SVM classification. The resultant training model was applied to identify pre-miRNA in the following study. In addition to the genome sequences, transcriptome sequences were also used to detect pre-miRNA sequences with the above identical method.

2.4 Prodiction of novel miRNAs using the ab initio approach

The overview workflow for detecting novel miRNAs is presented in Fig. 2. Briefly, all miRNAs, including the predicted miRNAs from the homology-based method, were first aligned to the large yellow croaker genome using BLASTN to locate their position with a word match size of 7 and an E-value cut-off of 2. Second, the miRNA cluster was defined with at least 3 miRNAs in a sequence window of 2 kb. To predict possible miRNA sequences in the range of the miRNA cluster, 80 and 40 bp were used as the pre-miRNA and step length, respectively. Pre-miRNA sequences conforming to the above seven criteria were taken to predict potential miRNA molecules.

2.5 GO and KEGG pathway annotation for putative miRNA target genes

3′ UTR sequences longer than 20 bp were extracted using a python script according to the gene annotation and used for target gene prediction. To reduce false positives, we combined the results of the prevalent miRNA target prediction tools of RNAhybrid and Miranda. The default score threshold and energy cut-off (score≥165 and energy≤-23) were used in miRanda to predict the miRNA target. The p value and E-value (P=0.1 and e=-23) were applied in RNAhybrid.

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used for the functional annotation of miRNA target genes. First, the miRNA targets genes were annotated by BLASTX against the NR database. Second, the best hit results by NR blast were used for GO and KEGG annotation with default settings in Blast2GO. We also used KAAS (KEGG Automatic Annotation Server) for KEGG annotation, and the results were combined with the output from Blast2GO.

2.6 Identification of miRNA-related SNP

Whole-genome resequencing data published in previous studies were used to call the SNP loci for large yellow croaker (Ao et al., 2015). The resulting SNP information was employed to analyse the SNP distribution around miRNA molecules. To investigate the characterization of the SNP distribution around miRNA, we estimated the density of SNPs in six regions (Gong et al., 2012), including miRNAs, pre-miRNAs, upstream and downstream 1 000 bp of pre-miRNA, and pre-miRNA upstream and downstream flanking regions (Fig. 3). We also randomly selected regions of 500 bp from the reference genome as the control region. The SNP density was defined as the number of SNPs per 1 000 bp.

2.7 The analysis of SNP-induced miRNA target gain and loss

The seed regions of miRNAs were more conserved than other regions of mature miRNAs because the SNP loci in miRNA target regions may influence the miRNA regulation cascades. To investigate how SNPs influence the miRNA function, miRNAs whose target regions contain SNP loci were considered in this study. As described above, two miRNA target tools, Miranda and RNAhybrid, were employed to predict the target genes for the wild-type and SNP-type miRNAs. For easy analysis of the effect induced by SNP in miRNAs, we defined the following four signs: WR (target genes of wild-type miRNAs processed by RNAhybrid), SR (target genes of SNP-type miRNAs processed by RNAhybrid), WM (target genes of wild-type miRNAs processed by Miranda), and SM (target genes of SNP-type miRNAs processed by Miranda). If one target exists in both WR and WM, but not in SM or SR, we defined this as target loss. If one target exists in both SM and SR, but not in WR or WM, we defined this as target gain. If one target exists in all four signs, we defined this as the target remain. All target genes for wild-type and SNP-type miRNA were classified into three groups, gained targets, lost targets and remained targets, for each miRNA.

3 RESULT AND DISCUSSION 3.1 Determination of the best E-value parameter for BLASTX

According to previous studies, the sequence features of the miRNA precursor, including the conformation of the stem loop secondary structure, low MFE, relative narrow range of the A+U content and MFEI, could be used to distinguish true miRNA molecules from other potential sequences; however, it is still necessary to eliminate coding sequences by BLASTX from the miRNA precursor sequences against the known protein database. In previous studies, an arbitrary E-value was chosen during miRNA classification. In our work, we attempted to estimate the optimized E-value for miRNA identification from receiver operating characteristic (ROC) analysis. To this end, we use 599 medaka miRNAs from NGS (Li et al., 2010) as the positive dataset and 436 predicted miRNA from medaka CDS sequences as the negative dataset. Both the positive dataset and negative datasets were evaluated BLASTX against the NR database. We combined two E-values of the BLASTX results as the input dataset of the ROC analysis. As shown in Fig. 4, the optimized E-value was -5.611, leading to a specificity of 0.609 and sensitivity of 0.733.

3.2 Identification of miRNAs in large yellow croaker according to a homology-based search

Most mature miRNAs are evolutionarily conserved among species. Animal miRNAs, including 1 275 species and 26 332 miRNAs, were selected from the latest miRBase. Large yellow croaker genome and transcriptome sequences were scanned to identify conservative miRNA genes using the prediction procedure, as demonstrated in Fig. 1. The whole miRNA detection processes included two steps. First, potentially conserved miRNAs were searched from genome and transcriptome sequences using BLASTN and preliminary filtration, resulting in 1 085 and 189 homologous miRNA from large yellow croaker genome and transcriptome sequences, respectively. Second, protein-coding and redundant miRNA were omitted from the potential miRNAs to reduce false positive results. After the filtration and SVM classification in the second step, the total numbers of potential miRNAs from the genome and transcriptome were reduced to 273 and 54 miRNAs, respectively. Redundant miRNAs were removed from total predicted miRNAs before we compared total predicted miRNAs with previous report miRNAs in large yellow croaker by using BLASTN, and 287 unique miRNAs were gained. The 287 identified miRNA sequences (Supplementary Table 2) covered 103 previously reported miRNAs (Qi et al., 2014; Huang et al., 2016) and 184 newly predicted miRNAs. According to the miRNA classification, all of the predicted miRNAs of large yellow croaker were grouped into 100 families (Supplementary Table 2). There were 184 newly miRNAs when compared with previous studies' results and 3 reasons can account for this. The first one is that the source data was different, in our study, we used the whole genome sequences and comprehensive transcriptome sequences while in previous research genome sequences or specific condition transcriptome sequences were used. The second one is that parameters may different for same software in different reports, for example, we use the ROC to find the best E-value for BLASTX. The last one which is the most important one is that different features were used to train SVM model when compared with previous researches.

3.3 Novel miRNA detection in the large yellow croaker genome using the ab initio approach

Previous reports revealed that miRNA genes were frequently grouped as clusters in chromosomes. For example, miRNA clusters containing more than three miRNA genes were found in D. melanogaster (Aravin et al., 2003); miRNA clusters containing let-7a-1 and let-7f-1 were identified in the human genome (Lagos-Quintana et al., 2001). Therefore, it is possible to discover novel miRNAs in the large yellow croaker genome around known miRNA clusters. Taking the miRNA predicted in this study and 112 miRNA genes previously reported (Qi et al., 2014) as known miRNA, 391 miRNA genes were left after removing redundant miRNA. To determine the miRNA gene locations, we mapped all known miRNAs in the large yellow croaker genome, and 384 (98.2%) of all known miRNA genes were successfully mapped to the best unique location.

To the best of our knowledge, the strict definition of a miRNA cluster is still lacking. Altuvia et al. utilized a fixed length with 3 kb to define a miRNA cluster (Altuvia et al., 2005), while other researchers adopted a more flexible length. For example, Chan et al. selected a range from 1 kb to 50 kb for a miRNA cluster (Chan et al., 2012). In this study, we chose 2 kb as the length of a miRNA cluster because a longer length may overestimate miRNA clusters in the large yellow croaker genome. A miRNA cluster was defined as the genomic region including at least two miRNA molecules. With the method described in the Methods section, 131 novel miRNAs were de novo predicted and are listed in Supplementary Table 3.

We combined homology results and ab initio results together and 418 miRNAs were gained. The number is very similar with previous report in Equus caballus (Zhou et al., 2009) with 407 miRNAs. Our results are more comprehensive when compared with the results in Equus caballu, as the fact we used both genome and transcriptome sequences as the input sequence. In the next future, new data on other species and the comparison with data obtained in closely related species will probably help to estimate if we can expect to find more miRNAs in the large yellow croaker. Besides, limits in our study were also exist and pri-miRNAs which only contain one miRNA precursor were considered while there are some pri-miRNAs which can produce more than one miRNA precursor (Bartel, 2004).

Predict miRNAs by this using computational method is efficient and low-cost, some work still need to be done in the future. In order to exam this workflow's accuracy, QPCR validations would be exerted. Furthermore, small RNA sequencing will be necessary in the future to confirm the biological significance of these predicted miRNAs as the fact that cost of sequencing technology is becoming cheaper and cheaper.

3.4 Features of the predicted miRNA genes

We calculated the features of the 418 predicted miRNAs and their precursor sequences from our research. In our study, the length of pre-miRNA varied from 56 to 187 nt, with an average of 90.96 nt, which was consistent with that of animal hairpin sequences downloaded from miRBase. The length of mature miRNAs ranged from 16 to 27 nt, with an average of 21.74 nt, showing a comparable length distribution with other animal mature miRNAs in the miRBase database (Supplementary Fig. 1).

We analysed the base composition at each position in mature miRNAs; the statistical results are shown in Fig. 5a. The nucleotide distributions at each position of large yellow croaker miRNA were in accordance with other animals (Fig. 5b). The base composition confirmed previous results that cytosine (C) is used less often than the other the three nucleotides at position 19, while adenine (A) is found at the smallest proportion at position 23.

miRNA precursors as well as other RNAs, mRNAs, tRNAs and rRNAs could fold into hairpin secondary structures. Therefore, the minimal free energy (MFE) and minimal free energy index (MFEI) are two important indicators to distinguish miRNAs from other small RNAs (Ng and Mishra, 2007). The MF E-value of the predicted pre-miRNAs ranges from -112.3 to -16.2 kcal/mol, with an average of -39.97 kcal/mol, and the average of MFEI was found to be 0.87 kcal/mol. The G+C content is another feature of miRNA precursors (Zhou et al., 2009) because a higher G+C content in miRNA precursors could influence the stability of the secondary structures. The G+C content of miRNA precursors varied from 30.36 to 68.97 with an average of 43.75.

3.5 Potential miRNA target gene prediction and functional annotation

miRNAs regulate gene expression by binding to the 3-UTR of target genes in various organisms. miRNA target gene identification is an important to better understand the role of miRNAs in cellular functions and gene regulation networks (Lai et al., 2016); however, false positives remain one of the most challenging issues during miRNA target prediction. We combined the 418 and 112 miRNAs, for a total of 530 miRNAs, from this and a previous study for further analysis. In this work, we predicted the potential targets of 530 miRNA genes by using two prevalent tools, miRanda and RNAhybrid. As a result, 16 142 and 23 510 miRNA target genes were predicted from miRanda and RNAhybrid, respectively. To reduce false positives, only target genes predicted by the two tools were used for the subsequent analysis. As a result, a total of 16 053 targets for 507 miRNA genes were predicted, with an average of 31 target genes for each miRNA. We only found that 23 miRNAs were not assigned to target genes in both tools, which might be attributed to the uncompleted reference genome and limitation of the target gene prediction tools.

The annotation of Gene Ontology (GO) uncovered the putative molecule functions of the target genes as well as the diverse biological processes that are regulated by related miRNAs. The 16 053 potential target genes were used for functional categorization and pathway analysis. The BLAST2GO suite was used to annotate the miRNA target genes by GO databases into three independent categories, cellular components, molecular functions, and biological processes. As a result, 14 609 target genes were successfully associated with GO terms. Furthermore, 13 284, 13 152 and 13 125 target genes were assigned to cellular components, biological processes and molecular functions, respectively. A more detailed annotation for each category is represented in Fig. 6. In the biological process category, the cellular processes (GO:0009987, 11357) and single-organism processes (GO:0044699, 10726) were the two most enriched terms. For the molecular function category, the majority consisted of binding (GO:0005488, 11052) and catalytic activity (GO:0003824, 5589). In the cellular component category, many miRNA target genes were involved in the cell part (GO:0044464, 11864), which was followed by cells (GO:0005623, 11915).

With the help of KEGG annotation, we can elucidate the functions of putative target genes related to primary metabolism and secondary metabolite biosynthetic pathways. KAAS (KEGG Automatic Annotation Server) was used for pathway annotation of genes with a BLAST output (see Methods for details). Then, 6 308 target genes were assigned to KO terms and classified into 341 pathways after mapping the target genes to the reference pathways. All pathways were divided into six categories according to the molecular interaction and reaction networks, including Metabolism, Genetic Information Processing, Environmental Information Processing, Cellular Processes, Organismal Systems and Human Diseases (Supplementary Fig. 2). With the whole genome genes as the background, we employed KEGG annotations for pathway enrichment. Those pathways represented the biological pathway influenced most by miRNA regulations (Table 1). The enriched pathways included many biological pathways that have been reported in previous studies, such as melanogenesis (Felicetti et al., 2008), dopaminergic synapse (Fiore et al., 2008) and the wnt signalling pathway (Gokhale et al., 2010).

3.6 SNP distribution around miRNA in the large yellow croaker genome

Animal miRNA genes exhibited high conservation among organisms; therefore, we speculated that mutations around miRNA re suppressed by natural selection. The SNP distributions around miRNA molecules were previously depicted in human and rice genomes (Saunders et al., 2007; Zhu et al., 2012; Liu et al., 2013); however, to the best of our knowledge, there is no research on SNPs around miRNA genes of large yellow croaker. To compare the SNP distribution among diverse regions in the large yellow croaker genome, SNP densities of six regions were calculated, including the miRNA, pre-miRNA (miRNA precursor) and pre-miRNA adjacent upstream and downstream 1 000 bp regions; pre-miRNA flanking regions and whole-genome random regions (Fig. 7). As we expected, the pre-miRNA flanking regions had the highest SNP density (~31.21 per kb), while the mature miRNA regions had the lowest SNP density (~10.41 per kb). The results were consistent with a previous investigation in the human genome, which also confirmed the reliability of miRNA detection in this study.

3.7 Target alteration caused by SNPs in miRNA seed regions

Previous investigations of polymorphisms in miRNA seed regions showed that miRNA SNPs might suppress gene expression by improving the binding affinity, while other mutations could lead to novel gene expression regulation. Alterations of the gene regulation cascade by polymorphisms in seed regions of miRNAs may significantly influence normal biological pathways in living cells; therefore, studies on SNPs in miRNA provided an important approach to probe the function of miRNAs in gene regulation pathways. For example, the T allele in miR-196a2 (rs11614913) could decrease the risk of breast cancer. By contrast, miR-499 (rs3746444) AG/GG could increase the risk of breast cancer (Hu et al., 2014). However, the biological influence analysis of polymorphisms in the seed regions of miRNAs has barely been reported for teleost.

miRNA target gene prediction mainly relied on miRNA seed region (nucleotides 2–8) complementation; therefore, SNPs in mature miRNA, especially the seed regions, may have a profound influence on the miRNA function and biogenesis. In this study, 22 SNPs were identified in the seed regions of 19 miRNAs. All miRNAs only had one SNP in the seed regions, except miR-7154-3p, miR-2025-5p and Pcr-mir-scaffold1370_34620. RNAhybrid and miRanda were employed to predict target genes for all mutated and wild-type miRNAs, except miR-449c-3p, which was too short for target gene prediction after mutation. All target genes for wild-type and SNP-type miRNA were classified into three groups for each miRNA, gained targets, lost targets and remaining targets. As shown in Fig. 8, we found that most of site mutations in putative miRNAs could greatly increase the target gene number. By comparing the remaining and gained target gene numbers for each miRNA, we found that nucleotide alterations of the miRNA seed regions can significantly enhance complementation with their target gene (Supplementary Table 4, t-test P-value=0.01). Meanwhile, the number of gained target genes was obviously higher than that of the remaining and lost targets, which was similar to the investigation in the rice genome, indicating that mutations in miRNA seed regions might dramatically influence the miRNA regulation pathway.

miRNA-related SNPs could impair or enhance miRNA processing and alter the sites of processing. Many scientists have studied miRNA-related SNPs especially in human diseases. For example, Liu et al. (2012) established a miRNA-related SNP database for human species; Hu et al. (2016) have found a SNP site associate with plasma triglycerides levels and coronary heart disease in miR-4217 target gene. However most of these researches are related about SNPs in miRNAs' target 3′UTR sequences or pre-miRNAs or pre-miRNAs. Studies associated with SNPs which located in miRNA seed regions are very rare. So it's hard to find a proper previous research to compare with our results. Although miR-4217 in Hu's findings play an important role in human heart disease, the SNPs' location is not the seed region. But our results can provide a potential research direction in the near future for other species.

4 CONCLUSION

Most computational investigations for miRNA identification rely on the sequence conservation of mature miRNA among species. However, the homology-based approach can only detect known miRNA molecules. In this study, we adopted an integrated strategy combining the homology-based and ab initio methods, utilizing the characterization of the miRNA cluster to predict miRNAs from the large yellow croaker genome and transcriptome. A total of 418 miRNAs, including 287 miRNAs from homology searching and 131 miRNAs from the ab initio approach, were identified. Then, 16 053 putative target genes were predicted for the identified miRNAs, which were further annotated by GO term and KEGG pathway analyses. The SNP distribution around miRNAs revealed that miRNA mutations were suppressed. Then, 22 SNP loci were identified in large yellow croaker miRNA seed regions, and we demonstrated that these mutations could dramatically influence miRNA-gene regulation cascades. This study not only provides abundant miRNA resources for large yellow croaker for subsequent investigations as well as offers a reference strategy for miRNA detection and analysis for other non-model organisms.

5 DATA AVAILABILITY STATEMENT

The genome sequence data and transcriptome sequencing data that support the findings of this study are available in NCBI (accession No. JRPU00000000) with the identifier https://doi.org/10.1371/journal.pgen.1005118 and NCBI (SRA accession No. of SRR1509885) with the identifier https://doi.org/10.1371/journal.pone.0124432.

6 COMPETING INTERESTS

The authors declare there are no competing interests.

Electronic supplementary material

Supplementary material (Supplementary Figs.1–2 and Tables 1–4) is available in the online version of this article at https://doi.org/10.1007/s00343-018-7053-5.