High-Throughput Sequencing of MicroRNA Transcriptome and Expression Assay in the Sturgeon, Acipenser schrenckii

Abstract

Sturgeons are considered as living fossils and have very high evolutionary, economical and conservation values. The multiploidy of sturgeon that has been caused by chromosome duplication may lead to the emergence of new microRNAs (miRNAs) involved in the ploidy and physiological processes. In the present study, we performed the first sturgeon miRNAs analysis by RNA-seq high-throughput sequencing combined with expression assay of microarray and real-time PCR, and aimed to discover the sturgeon-specific miRNAs, confirm the expressed pattern of miRNAs and illustrate the potential role of miRNAs-targets on sturgeon biological processes. A total of 103 miRNAs were identified, including 58 miRNAs with strongly detected signals (signal >500 and P≤0.01), which were detected by microarray. Real-time PCR assay supported the expression pattern obtained by microarray. Moreover, co-expression of 21 miRNAs in all five tissues and tissue-specific expression of 16 miRNAs implied the crucial and particular function of them in sturgeon physiological processes. Target gene prediction, especially the enriched functional gene groups (369 GO terms) and pathways (37 KEGG) regulated by 58 miRNAs (P<0.05), illustrated the interaction of miRNAs and putative mRNAs, and also the potential mechanism involved in these biological processes. Our new findings of sturgeon miRNAs expand the public database of transcriptome information for this species, contribute to our understanding of sturgeon biology, and also provide invaluable data that may be applied in sturgeon breeding.

About the Authors:

Lihong Yuan

Contributed equally to this work with: Lihong Yuan, Xiujuan Zhang

Affiliation: Guangdong Public Laboratory of Wild Animal Conservation and Utilization, Guangdong Entomological Institute/South China Institute of Endangered Animals, Guangzhou, China

Xiujuan Zhang

Contributed equally to this work with: Lihong Yuan, Xiujuan Zhang

Affiliation: Guangdong Public Laboratory of Wild Animal Conservation and Utilization, Guangdong Entomological Institute/South China Institute of Endangered Animals, Guangzhou, China

Linmiao Li

Affiliation: Guangdong Public Laboratory of Wild Animal Conservation and Utilization, Guangdong Entomological Institute/South China Institute of Endangered Animals, Guangzhou, China

Haiying Jiang

Affiliation: Guangdong Public Laboratory of Wild Animal Conservation and Utilization, Guangdong Entomological Institute/South China Institute of Endangered Animals, Guangzhou, China

Jinping Chen

* E-mail: [email protected]

Affiliation: Guangdong Public Laboratory of Wild Animal Conservation and Utilization, Guangdong Entomological Institute/South China Institute of Endangered Animals, Guangzhou, China

Introduction

Sturgeons (order: Acipenseriformes, infraclass: Chondrostei), are believed to have separated from other teleosts over 250 MYA, and are referred as living fossils providing a key phylogenetic position for evolutionary studies on vertebrates [1]. As different degrees of ploidy are resulted from the multiple and independent duplication events [2], [3], sturgeon species have been divided into two groups (approximately 120 and 240 chromosomes), and the division of the species into two groups of either diploid/tetraploid or tetra/octoploid is still unclear [4], [5]. A previous study showed that microRNAs (miRNAs) in Arabidopsis thaliana evolved with a process of genome-wide duplication followed by dispersal and diversification, and the duplicated copies of miRNAs may acquire new functionality during their evolution [6]. A later study proved that almost all of miRNAs of Rainbow trout (Oncorhynchus mykiss) have been retained as duplicated copies [7]. These results indicated that, by chromosome duplication, new miRNAs in sturgeon may be evolved, and thus reflect new functions. Hence, in-depth study of sturgeon miRNAs would improve the knowledge of sturgeon ploidy, and also evolutionary history of sturgeons and vertebrates.

MiRNAs are small RNA molecules (18–25 nt), and have been identified to play key roles in directing transcriptional and post-transcriptional activity of mRNAs, and are thus involved in the regulation of multiple biological processes such as differentiation and development, immune response, reproductive system development and gametogenesis [8], [9]. The study of the sturgeon miRNAs and their interactions with target mRNAs will provide further insight into these physiological processes of sturgeon. RNA-seq, an ultrahigh-throughput sequencing technique has greatly increased our understanding on the complexity of eukaryotic mRNA and small RNA transcriptomes [10], including those of non-model species, which have neither genomic background nor miRNA data in miRbase [11], [12].

In the present study, we carried out the first high-throughput sequencing analysis on miRNAs of five key tissues (liver, spleen, muscle, heart and brain) of the Amur sturgeon Acipenser schrenckii, an endangered and important economic sturgeon species, using the RNA-seq technique on the Illumina TruSeq sequencing platform. Combining with the expression level validation of miRNAs by microarray and stem-loop real-time PCR, our study aims to discover sturgeon-specific miRNAs, investigate the expression pattern and illustrate the potential role of miRNAs and their targets on sturgeon biological processes.

Materials and Methods

Ethics statement

The protocol was approved by the Committee on the Ethics of Animal Experiments of the Guangdong Entomological Institute, which also incorporates the South China Institute of Endangered Animals. Sturgeon individuals were immerged in the water with 10⁻⁴ (v/v) Eugenol about 1–3 minutes for euthanasia, following the AVMA guidelines (2013) for use [13]. All efforts were made to minimize suffering.

Sample and RNA preparation

The five tissues (liver, spleen, muscle, heart, brain) of a 5-month-old Amur sturgeon, Acipenser schrenckii, from the Engineering and Technology Center of Sturgeon Breeding and Cultivation of Chinese Academy of Fishery Science (Beijing, China) were collected. Total RNA were extracted from five tissue samples separately with RNAiso reagent (TaKaRa, Japan) according to the manufacturer’s instructions. RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies), and RNA purity was assessed using the Nano Photometer spectrophotometer (IMPLEN). Total RNA of five tissues (3 ug each) were pooled, and RNA integrity was inspected using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies). The pooled RNA sample with RNA Integrity Number (RIN) = 8.3 met the needs of TruSeq transcriptome/small RNA library construction and sequencing, which hereafter referred to as ASY transcriptome/small RNA library.

Transcriptome references sequencing

As no genomic sequences specific to A. schrenckii are available on public database, we firstly carried out the de novo transcriptome sequencing and assembling by Illumina TruSeq platform. 3 ug pooled RNA was used as mRNA library construction using Illumina TruSeq RNA Sample Preparation Kit (Illumina) following manufacturer’s recommendations. Briefly, mRNA was purified from 3 ug pooled RNA by using poly-T oligo-attached magnetic beads. After the first and second strand cDNA synthesizing, DNA fragments were converted into blunt ends, adenylated the 3′ ends, and then ligated with Illumina PE adapter oligonucleotides for hybridization. Then, cDNA fragments with length >200 bp were purified with AMPure XP system (Beckman), and those ligated with adapters on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 10 cycles PCR reaction, and the products were purified again by AMPure XP system and quantified by Agilent 2100 bioanalyzer. Subsequently, the cluster of index-coded samples was generated using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) and sequenced on an Illumina Hiseq 2000 platform. Finally, 100 bp paired-end reads were generated. After removing the reads with adapters, any reads containing ‘n’ (>10%), low quality reads (sQ≤5) and the redundant reads, the remaining clean reads were assembled by TRINITY method [14], and then the redundant contigs were screened by CAP3 [15]. Finally, the unigenes were searched against Nr database (NCBI non-redundant protein sequences) by Blast2GO [16], and the orthologs were used as the reference sequences. All cDNA data series were submitted to NCBI Sequence Read Archive (SRA) database with accession number SRR1131121.

Construction and high-throughput sequencing of Small RNA library

According to the protocol of Illumina TruSeq Small RNA Sample Preparation Kit (Illumina), 3 ug pooled RNA was used as small RNA library construction. In brief, RNA bands around 20–30 bp were separated and purified by 6% TBE PAGE gel and subsequently bound to 3′ and 5′ end adapters in two separated subsequent steps, which followed by PAGE gel purification. After the first strand cDNA synthesizing by random oligonucleotides and SuperScript II and amplifying by PCR, DNA fragments ligated with adapters on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 12 cycles PCR reaction, and the products of 145 bp to 160 bp (with adaptors on both sides) were separated by PAGE gel, and quantified by Agilent 2100 bioanalyzer. Then, the cluster of index-coded samples was generated using TruSeq SE Cluster Kit v3-cBot-HS (Illumina) and sequenced on an Illumina Hiseq 2000 platform. Finally, 50 bp single-end reads were generated. All small RNA data series were submitted to SRA database with accession number SRR1129970.

Filter of small RNA reads and microRNAs identification

After removing the unclean reads (the adapters, low quality reads, reads containing ‘n’, and redundant reads), clean unique reads were mapped onto the A. schrenckii transcriptome reference sequences using the program Bowtie [17] with no mismatch. Perfectly mapped reads were scanned against the Metazoa mature microRNA (miRNA) of Sanger miRBase (Release 19) [18] to identify the orthologs of known miRNAs. Then, the non-conserved unique reads were screened against Rfam (http://rfam.sanger.ac.uk/) [19] and RepeatMasker (http://www.repeatmasker.org/) [20] successively using the program Bowtie to filter the sequences originating from rRNA, tRNA, snRNA, snoRNA and repetitive elements.

The potential miRNA reads, which were unannotated small RNA tags and could be mapped onto the transcriptome reference sequences, were analyzed by miREvo [21] and mirdeep2 [22] for the prediction of Dicer cleavage site, the assay of secondary structure and the minimum free energy. Finally, the potential miRNA candidates were submitted to miRBase again, and the precursors (hairpins) of potential miRNAs that passed MirCheck [23] were manually inspected the canonical structure of miRNAs in order to remove the false prediction. The base bias of mature miRNAs on the first nucleotide position with certain length and on each position of all identified miRNAs (>100 reads) were calculated, respectively.

MiRNA microarray and data analysis

We used another 5-month-old Amur sturgeon individual from the Engineering and Technology Center of Sturgeon Breeding and Cultivation of Chinese Academy of Fishery Science (Beijing, China) to validate the expression of 103 miRNAs identified by Illumina TruSeq sequencing. The five tissues (liver, spleen, muscle, heart, brain) were sampled and total RNA were extracted as described above. 4 ug total RNA of each sample was used to hybridize with microarray chip.

MiRNA microarray was manufactured by LC Sciences (China), and each miRNA probe has five replicates. The chip was hybridized with RNAs of five Amur sturgeon tissues, which were 3′-extended with a poly (A) tail and ligated with an oligonucleotide tag. Hybridization was performed overnight on a µParaflo microfluidic chip using a micro-circulation pump (Atactic Technologies) [24]. Single-color labeling (Cy3) and hybridization of total RNA were performed according to the manufacture’s protocol with no modification. Microarray results were extracted using a laser scanner (GenePix 4000B, Molecular Device) and digitized using Array-Pro image analysis software (Media Cybernetics). Raw data were subtracted by the background matrix, and then normalized using a LOWESS (Locally-weighted Regression) method to remove system related variations, including sample amount variations and signal gain differences of scanners, and thus faithfully reveal the biological variations [25]. The transcripts will be defined as detectable when their signal intensity higher than 3× (background standard deviation), spot CV [(standard deviation)/(signal intensity)]<0.5, and 50% of the repeating probes are meet the first two criteria. Finally, the expression of miRNAs was determined with the criterion of signal intensity and P value. The microarray data were deposited in a database (ArrayExpress, GEO) with accession number GSE57102.

Quantitative miRNA real-time PCR assay

Stem-loop reverse transcription (RT) real-time PCR was used to quantify the expression of ten mature miRNAs in five tissues (brain, heart, liver, spleen and muscle) of three 5-month-old Amur sturgeon individuals. U6 snRNA, which has relatively stable expression in most of the tissues, was used as the endogenous control [26]. Briefly, 500 ng of total RNA of each sample was reverse-transcribed with miRNA-specific stem-loop RT primers using the First-strand cDNA Synthesis Kit (Thermo Scientific Fermentas). The reactions were incubated at 42°C for 60 min, at 70°C for 15 min and then held at 4°C. Real-time PCR was performed in triplicate wells using Strategene Mx3000P (Agilent Technologies company, American) according to the protocol. In a 20 µl reaction mixture, 2.0 µl of cDNA was used as template, with 10 µl of SYBR Select Master Mix (Applied Biosystems, Carlsbad, USA), 1.0 µl of specific forward primer, and 1.0 µl of universal primer, with the following program: 50°C for 2 min for UDG (Heated-labile Uracil-DNA Glycocasylase) activation and then 95°C for 2 min, followed by 40 cycles of 95°C for 15 s, 60°C for 30 s and 72°C for 30 s. Primers used in this study were listed in S1 Table. The 2^−ΔCT method was used to calculate the relative expression (versus U6 snRNA).

MiRNA targets prediction and annotation

Hereafter, we only considered the miRNAs detected by microarray with the criteria of Signal >500 and P≤0.01 for further analysis.

For miRNAs target gene prediction, we first predicted the Open Reading Frame (ORF) of Amur sturgeon reference sequences, identified the orthologous mRNAs and then predicted the 3′UTR by searching against the vertebrate genomic database in GENSCAN (http://genes.mit.edu/GENSCAN.html). The 3′ UTR sequences of the orthologs were trimmed and analyzed by using the Perl scripts of both TargetScan v6.2 with context score percentile ≥50 (http://www.targetscan.org/) and miRanda v3.3a with Max_Energy ≤ −20 (http://www.microrna.org/microrna/home.do) for the target gene prediction. Then, the target genes were annotated by mapping to Gene Ontology (GO) database (http://www.geneontology.org/) and KEGG pathways database (http://www.genome.jp/kegg/) by BLASTX at E values ≤1e-5 [27]. Finally, the enriched functional groups or pathways among miRNA putative targets were identified with P<0.05.

Results

Transcriptome references sequencing

A total of 5.89×10⁷ reads were sequenced from the ASY transcriptome library with error rate of 0.04%, Q30 of 87.52% and GC content of 49.61%. Total 5.01×10⁷ (85.07%) clean reads remained after removing the low quality and contaminant reads (S1 Figure). After reads assembly, removing the redundancy and annotation of unique sequences, a total of 148,817 unigenes (N₅₀ = 1599) specific to Acipenser schrenckii were obtained, including 41,378 protein-coding sequences.

Sequencing and statistics of small RNA reads

A total of 1.67×10⁷ reads were sequenced from the ASY small RNA library with error rate of 0.01%, Q30 of 94.14% and GC content of 53.61%. Then a total 1.35×10⁷ high-quality small RNA reads were obtained after removing the ambiguous reads (Table 1). The size distribution and frequency percentage of small RNA reads are shown in S2 Figure, and in them, the potential miRNA reads (21–24 bp) were the major part (about 44%).

[Figure omitted. See PDF.]

Table 1. Classification of small RNA reads of Acipenser schrenckii.

https://doi.org/10.1371/journal.pone.0115251.t001

After mapping to the A. schrenckii transcriptome reference sequences, 8.37×10⁶ perfectly matched small RNA reads remained. A total of 5.5×10⁵ reads, were identified by searching against miRBase, corresponding to 6.58% perfectly matched small RNA reads (Table 1). Subsequent small RNAs filter showed that other non-coding RNAs (rRNA, tRNA, snRNA and snoRNA), repeat sequences and unknown genomic regions were about 21.73%, 0.16% and 70.78%, respectively. In those, the repeat sequences were further inspected, and DNA (DNA transposons), minisatellite, LTR (Long Terminal Repeats) and LINE (Long INterspersed Elements) are the most abundant parts (S3 Figure). Finally, 6.4×10⁴ (0.76%) potential novel miRNA reads specific to A. schrenckii were detected by miREvo (Table 1).

MiRNA identification

After the secondary structures analysis and manual examination, we identified 75 miRNA precursors with the stem-loop hairpin structures characteristic to miRNAs precursors, corresponding to 52 unique mature miRNAs named with ASY-X where Xs are numeric numbers (Table 2 and S2 Table). To predict putative miRNAs, we applied miRDeep2 program, which incorporates the position and frequency of small RNAs with the secondary structure of miRNA precursor and can discover novel miRNAs [28], [29], on the potential novel miRNA reads specific to A. schrenckii. We obtained 51 unique miRNA precursors, corresponding to 51 uniquely putative novel miRNAs named ASY-novel-X (Table 2 and S2 Table). In those, three putative novel miRNAs were found in miRBase (ASY-novel-5 and ASY-novel-51 are identical with pol-miR-144-5p and pol-miR-133-5p respectively; ASY-novel-37 is similar to dre-miR-1306 with one mismatch). Thus, these three miRNAs were classified into the conserved miRNAs and named as ASY-miR-144-5p, ASY-miR-133-5p and ASY-miR-1306. The secondary structure of most abundant putative novel miRNA (ASY-novel-46) is shown in Fig. 1.

[Figure omitted. See PDF.]

Figure 1. Secondary structure of the novel miRNA (ASY-novel-46) specific to Acipenser schrenckii as predicted by miRDeep2.

https://doi.org/10.1371/journal.pone.0115251.g001

[Figure omitted. See PDF.]

Table 2. Top 20 highest expressed miRNAs of Acipenser schrenckii identified by Illumina sequencing.

https://doi.org/10.1371/journal.pone.0115251.t002

The base-bias analysis indicated the strong preference of nucleotide utility in both conserved and putative novel miRNAs (Fig. 2A and 2B). Among conserved miRNAs, the preference of first base utility shows a close correlation with length of miRNAs. The miRNAs with the length of 18 bp prefer the nucleotide C at the first base position, those of 19–24 bp favor the U (Fig. 2A, top). Similarly, the first base utility of putative novel miRNAs prefer the A in those with 18 bp length and the U in those with the length of 20 bp, whereas others with length of 21–23 bp favor the use of C/U or A/U (Fig. 2A, bottom). Moreover, we found that the nucleotide utility shows difference on each site position between conserved and putative novel miRNAs. For conserved miRNAs, the nucleotide U has the most frequency among sites, then followed by G, A and C (Fig. 2B, top). Whereas, the nucleotide A has the highest utility rate among putative novel miRNAs, and other three nucleotides (G, C and U) share the relative equal utility (Fig. 2B, bottom).

[Figure omitted. See PDF.]

Figure 2. Base bias analysis of the first nucleotide of Acipenser schrenckii miRNAs.

(A) Base bias percentage on the first position of conserved miRNAs (top) and novel miRNAs specific to A. schrenckii (bottom) sized from 18 ∼26 bp. The x-axis indicates the length of miRNAs. The y-axis indicates the percentage of four nucleotide acids. The numbers on the columns are the total numbers of miRNAs with specific length of each nucleotide. (B) Base bias percentage from the first to the 26th base pairs of conserved miRNAs (top) and novel miRNAs (bottom). The x-axis indicates the location of base pairs. The y-axis indicates the percentage of four nucleotide acids.

https://doi.org/10.1371/journal.pone.0115251.g002

Expression pattern assay of miRNAs by microarray and real-time PCR

We used the independent microarray platform to validate the expression level of 103 miRNAs obtained by Illumina TruSeq sequencing. A total of 87 miRNAs were detected by microarray in at least one of five tissues (S3 Table), and in those, 57 miRNAs have strong detected signal with the criterion of Signal >500 and P<0.01; 14 miRNAs had lower signal with P<0.01, and Signal <500 (S3 Table). Moreover, we found ASY-miR-21, the highest expressed miRNA in TruSeq, had detectable signal with Signal = 2.5×10⁵ and P = 0.0108; ASY-novel-1 with high expressed level in TruSeq had low detectable signal (Signal <500), and P<0.01 (S3 Table).

Further analysis indicated that, in 58 miRNAs with high detected signal (Signal >500 and P≤0.01), 21 miRNAs showed co-expression in all five tissues, and 16 miRNAs (one in liver, three in spleen, three in muscle and nine in brain) highly expressed in one specific tissue (Table 3). The expressed pattern of 58 miRNAs (Signal >500 and P≤0.01) in five tissues indicated that miRNAs were mainly separated into two clades, clade I, 46 of these are expressed at low levels in most of tissues, however, some of these are highly expressed in specific tissues such as brain, muscle and heart; and clade II, 12 miRNAs are very highly expressed in most of five tissues (Fig. 3).

[Figure omitted. See PDF.]

Figure 3. Expressed pattern of miRNAs in five tissues of Acipenser schrenckii.

The expression of 58 miRNAs detected by microarray (Signal >500 and P≤0.01) are reflected by Log-normalized intensities. Heat map represents the miRNAs which were clustered into two clades based on their expression in tissues.

https://doi.org/10.1371/journal.pone.0115251.g003

[Figure omitted. See PDF.]

Table 3. Details of miRNAs co-expressed in five tissues and specifically high-expressed in one tissue obtained by microarray.

https://doi.org/10.1371/journal.pone.0115251.t003

We randomly tested the expression of ten miRNAs (seven conserved miRNAs and three novels specific to A. schrenckii) in five tissues of three sturgeon individuals respectively by stem-loop real-time PCR assay. All these miRNAs have high expressed levels detected by both Illumina TrueSeq sequencing and microarray, with the exception of ASY-novel-1 which shows low expression level in microarray. Real-time PCR results supported the expression pattern of miRNAs obtained by microarray, with the exception of ASY-miR-101a, which showed much lower expression level in heart than that seen in the microarray (Fig. 3 and Fig. 4).

[Figure omitted. See PDF.]

Figure 4. Profiling of ten sturgeon miRNAs expression by real-time PCR.

Figure shows the relative expression of miRNAs (vs. U6 snRNA) and transformed by Log2. Data are presented as mean±SD (standard deviation) and n = 3.

https://doi.org/10.1371/journal.pone.0115251.g004

Target prediction and annotation

To better understand the functions of miRNAs, the putative targets of 87 miRNAs, which detected by both Illumina TruSeq sequencing and microarray were predicted by TargetScan and miRanda based on the A. schrenckii transcriptome. Firstly, in 148,817 of A. schrenckii unigenes, a total of 68435 genes, orthologous to vertebrate mRNAs and with the best hit, were remained. In these, 14265 of 68435 A. schrenckii unigenes were identified as the targets of miRNAs by both TargetScan (context score percentile≥50) and miRanda (Max_Energy ≤ −20). A total of 3391 GO terms and 228 KEGG pathways regulated by 87 miRNAs were also identified. In those, 369 enriched functional groups and 37 enriched functional KEGG pathways, which were mainly regulated by 58 high-expressed miRNAs, were obtained with P<0.05 (S4 Table). Most of the enriched functional pathways or groups were involved in cell development, signal transduction, metabolic and immune processes.

Discussion

Our study shows that total 58 of 103 microRNAs (miRNAs) were confirmed by both Illumina TruSeq sequencing and microarray in Amur sturgeon Acipenser schrenckii, including 21 miRNAs co-expressed in all five tissues and 16 which have tissue-specific expression. According to the functional annotation and the enrichment analysis of putative targets, these miRNAs are mainly involved in development, metabolism, immune response and gametogenesis.

For small RNA filtration and miRNA annotation, we firstly performed the transcriptome reference sequencing by Illumian TruSeq platform. Total 148,817 unigenes, which are assembled by 5.01×10⁷ clean reads and specific to A. schrenckii, were obtained. By small RNA sequencing, we obtained a total of 1.35×10⁷ high-quality small RNA reads and of them, about 44% are the potential miRNA reads with 21–24 bp in length (Table 1 and S2 Figure). After mapping to the A. schrenckii transcriptome reference sequences, about 7% reads were identified as potential miRNAs by searching against miRBase and miREvo detection, and the unknown genomic regions were the major part about 70%. The small mapping proportion of potential miRNA reads and a large number of unknown genomic sequences in small RNAs reads of sturgeon are not unique, whereas have been observed in some organisms, such as ∼3% vs. ∼94% in sea cucumber [12] and ∼1% vs. ∼97% in swithgrass [30]. This result may because of the lack of genomic background and limited transcriptome information specific to the Amur sturgeon. Moreover, another reason there is a low amount of small RNA reads is the degradation of RNA sample, which is indicated by the high proportion of rRNA (∼21%). In this study, we conducted the TruSeq sequencing with RNA pool of five tissues (RIN = 8.3), whereas, we did not test the RIN of five RNA samples separately before pooling samples. Thus, there is a chance that some sample degradation may have happened before pooling.

Our data also revealed a large number of sturgeon specific miRNAs, and provided candidates for further study of sturgeon biology. A total of 103 mature miRNAs were identified by TruSeq, including 55 conserved miRNAs and 48 novel miRNAs (S2 Table). Further analysis suggested the strong preference of nucleotide utility of sturgeon miRNAs, and also the great difference of nucleotide utility between conserved and novel miRNAs, especially in the first base position (Fig. 2). Previous studies have shown that the first nucleotide at the 5′ end of miRNA is considered key for strand selectivity of Dicer-mediated cleavage [31]. The strong preference of nucleotide utility on the first base position of sturgeon miRNAs may affect the Dicer cleavage, and thus, the target recognition [31], [32]. Further comparing the preference of nucleotide utility of sturgeon miRNAs with those of other organisms will increase the understanding on the interaction of sturgeon miRNAs-mRNAs.

Of 103 sturgeon miRNAs, 87 were detected by microarray validation including 58 miRNAs with strongly detected signals (S3 Table). Moreover, of those miRNAs with strongly detected signals, we randomly tested the expression pattern of 10 miRNAs in five sturgeon tissues by real-time PCR, and obtained consistent results with that of microarray (Fig. 3 and Fig. 4). In addition, we found that 21 miRNAs were co-expressed in all five tissues, and this suggested the crucial role of them for sturgeon physiological processes (Table 3). In them, 8 miRNAs (miR-21, miR-30c, miR-126-3p, let-7c, let-7e, miR-128, miR-20b and miR-181b) had been shown to be related with gametogenesis [8], [33]. Furthermore, 16 miRNAs were tissue-specifically expressed (Table 3), indicated the particular roles of these miRNAs in the related tissues. MiR-144, with the co-regulation of miR-451, has been proven to be involved in many diseases, such as anemia severity of sickle cell disease [34], cancer [35], brain aging and spinocerebellular ataxia pathogenesis [36]. The high-expression of miR-144 specifically in spleen, an immune and hematopoietic organ, and also miR-451 in all five tissues of sturgeon suggest the key role of miR-144/451 in sturgeon immune system. The specific high-expression of miR-133a in muscle, which has function in the proliferation and differentiation of cardiomyocytes, bronchial smooth muscle and related diseases [37]–[39], implies the role of miR-133a in sturgeon physiological processes. Moreover, studies showed that miR-9 plays role in neural diseases and cancers [40]–[42], miR-219 regulates neural precursor maintenance and specification [43], miR-34a is a tumor suppressor [44], miR-129 promotes apoptosis [45], and miR-7a regulates the neuronal excitability [46], and all these miRNAs were found to be specifically expressed in sturgeon brain. The identification of tissue-specifically expressed miRNAs, combining with the co-expressed miRNAs in all five tissues, provides clues to further study the molecular mechanism of sturgeon physiological processes.

KEGG pathway and GO annotation analyses could provide a better understanding of the potential functions of miRNAs by illustrating the function of target mRNAs. Because of the absence of the genomic information and incomplete transcriptome annotation of sturgeon, we used the genomic sequences of all vertebrate species deposited in the online gene analysis software (GENSCAN) as reference to predict the target gene of Amur sturgeon miRNAs. The GO terms and KEGG pathways identified here (S4 Table), especially the enriched functional groups and pathways which regulated by 58 miRNAs with strongly detected signal, provide us with guidance to purposely comb the miRNAs and putative mRNAs from the complex gene database networks, and can be used in the future for sturgeon aquaculture.

With the progress of the sturgeon genomic and transcriptome information, or for other closer species that become available, the accuracy of sturgeon gene annotation will be increased and much more detailed information of A. schrenckii will be uncovered from our mRNA and miRNA transcriptome datasets. Furthermore, compared with other teleost aquaculture species, the long lifespan and time to breeding capability of the sturgeon results in high maintenance costs (eg. space, energy consumption and water purification). MiRNAs that may be involved in cell differentiation and development, signal transduction, gametogenesis, metabolic and immune processes, and their targets provide opportunity for early sex determination and also to select individuals with rapid growth and disease resistance for breeding purposes.

Conclusions

This study reveals the first sturgeon miRNAs profile, and the findings advance the understanding of sturgeon biology and are valuable for sturgeon fishery and conservation. However, the validation of the relationship between sturgeon miRNAs and target mRNAs in the regulation of specific physiological processes needs further biologically experimental evidences.

Supporting Information

[Figure omitted. See PDF.]

S1 Figure.

Overview of Acipenser schrenckii transcriptome sequencing reads.

https://doi.org/10.1371/journal.pone.0115251.s001

(PNG)

S2 Figure.

The sequence length distribution and frequence percentage of small RNA reads of Acipenser schrenckii. The x-axis indicates the length of small RNA reads. The y-axis indicates the percentage of small RNA reads with specific length. Different color suggests different type of small RNAs.

https://doi.org/10.1371/journal.pone.0115251.s002

(PNG)

S3 Figure.

Classification of repeat sequences of Acipenser schrenckii small RNA library. Ambi: ambiguous reads; RC, rolling circle; LINE, Long INterspersed Elements; SINE, Short INterspersed Elements; LTR, Transposable elements with Long Terminal Repeats; DNA, DNA transposons. +, sense strand; -, anti-sense strand.

https://doi.org/10.1371/journal.pone.0115251.s003

(PNG)

S1 Table.

Forward, stem-loop and universal primers used to amplify miRNAs and U6 snRNA in real-time PCR.

https://doi.org/10.1371/journal.pone.0115251.s004

(DOC)

S2 Table.

Details of mature miRNAs and hairpins screening by Illumina. (A) 55 conserved mature miRNAs; (B) 75 conserved miRNA hairpins; (C) 48 mature novel miRNAs specific to A. schrenckii; (D) 51 noval miRNA hairpins.

https://doi.org/10.1371/journal.pone.0115251.s005

(XLSX)

S3 Table.

Summary of 87 miRNAs identified by both Illumina TruSeq sequencing and microarray. RPM, Reads Per Million.

https://doi.org/10.1371/journal.pone.0115251.s006

(XLS)

S4 Table.

Details of GO annotation and KEGG pathway analysis. The enriched functional groups and pathways were identified by Fisher’s Exact Test (P<0.05). S gene number: Number of Significant genes matched to single GO term or KEGG pathway; TS gene number: Total number of Significant genes matched to GO terms or KEGG pathways; B gene number: Number of genes matched to single GO term or KEGG pathway; TB gene number: Total number of genes matched to GO terms or KEGG pathways.

https://doi.org/10.1371/journal.pone.0115251.s007

(XLS)

Acknowledgments

We thank Mr. Wenhua Wu and Dr. Ying Zhang for help with sturgeon sampling, and Bronwyn McAllan and two anonymous reviewers for commenting on the manuscript.

Author Contributions

WH.W and Y.Z help for sample collection. Conceived and designed the experiments: LHY JPC. Performed the experiments: LHY XJZ LML HYJ. Analyzed the data: LHY XJZ. Contributed reagents/materials/analysis tools: LHY JPC. Contributed to the writing of the manuscript: LHY JPC

Citation: Yuan L, Zhang X, Li L, Jiang H, Chen J (2014) High-Throughput Sequencing of MicroRNA Transcriptome and Expression Assay in the Sturgeon, Acipenser schrenckii. PLoS ONE 9(12): e115251. https://doi.org/10.1371/journal.pone.0115251

References

1. Bemis WE, Findeis EK, Grande L (1997) An overview of Acipenseriformes. Sturgeon Biodiversity and Conservation 17:5–71.

2. Fontana F, Congiu L, Mudrak VA, Quattro JM, Smith TI, et al. (2008) Evidence of hexaploid karyotype in shortnose sturgeon. Genome 51:113–119.

3. Ludwig A, Belfiore NM, Pitra C, Svirsky V, Jenneckens I (2001) Genome duplication events and functional reduction of ploidy levels in sturgeon (Acipenser, Huso and Scaphirhynchus). Genetics 158:1203–1215.

4. Fontana F, Lanfredi M, Kirschbaum F, Garrido-Ramos MA, Robles F, et al. (2008) Comparison of karyotypes of Acipenser oxyrinchus and A. sturio by chromosome banding and fluorescent in situ hybridization. Genetica 132:281–286.

5. Birstein VJ, Vasiliev VP (1987) Tetraploid-octoploid relationships and karyological evolution in the order Acipenseriformes (Pisces) karyotypes, nucleoli, and nucleolus-organizer regions in four acipenserid species. Genetica 72:3–12.

6. Maher C, Stein L, Ware D (2006) Evolution of Arabidopsis microRNA families through duplication events. Genome Research 16:510–519.

7. Berthelot C, Brunet F, Chalopin D, Juanchich A, Bernard M, et al. (2014) The rainbow trout genome provides novel insights into evolution after whole-genome duplication in vertebrates. Nature Communications 5:3657.

8. Niu Z, Goodyear SM, Rao S, Wu X, Tobias JW, et al. (2011) MicroRNA-21 regulates the self-renewal of mouse spermatogonial stem cells. Proceedings of the National Academy of Sciences 108:12740–12745.

9. Bartel DP (2004) MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116:281–297.

10. Wang Z, Gerstein M, Snyder M (2009) RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10:57–63.

11. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, et al. (2013) De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nature Protocols 8:1494–1512.

12. Chen M, Zhang X, Liu J, Storey KB (2013) High-throughput sequencing reveals differential expression of miRNAs in intestine from sea cucumber during aestivation. PLoS One 8:e76120.

13. Leary S, Underwood W, Anthony R, Cartner S, Corey D, et al. (2013)AVMA guidelines for the euthanasia of animalsedition: 2013 edition.

14. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, et al. (2011) Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology 29:644–652.

15. Huang X, Madan A (1999) CAP3: A DNA sequence assembly program. Genome Research 9:868–877.

16. Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, et al. (2005) Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21:3674–3676.

17. Langmead B, Trapnell C, Pop M, Salzberg SL (2009) Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology 10:R25.

18. Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ (2008) miRBase: tools for microRNA genomics. Nucleic Acids Research 36:D154–158.

19. Burge SW, Daub J, Eberhardt R, Tate J, Barquist L, et al. (2013) Rfam 11.0: 10 years of RNA families. Nucleic Acids Research 41:D226–232.

20. Tarailo-Graovac M, Chen N (2009) Using RepeatMasker to identify repetitive elements in genomic sequences. Current Protocols in Bioinformatics25: 4.10. 1–14.10:14.

21. Wen M, Shen Y, Shi S, Tang T (2012) miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics 13:140.

22. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N (2011) miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Research 40:37–52.

23. Jones-Rhoades MW, Bartel DP (2004) Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell 14:787–799.

24. Gao X, Gulari E, Zhou X (2004) In situ synthesis of oligonucleotide microarrays. Biopolymers 73:579–596.

25. Bolstad BM, Irizarry RA, Astrand M, Speed TP (2003) A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19:185–193.

26. Peltier HJ, Latham GJ (2008) Normalization of microRNA expression levels in quantitative RT-PCR assays: identification of suitable reference RNA targets in normal and cancerous human solid tissues. Rna 14:844–852.

27. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, et al. (1997) Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 25:3389–3402.

28. Friedländer MR, Chen W, Adamidi C, Maaskola J, Einspanier R, et al. (2008) Discovering microRNAs from deep sequencing data using miRDeep. Nature Biotechnology 26:407–415.

29. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N (2012) miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Research 40:37–52.

30. Xie F, Stewart CN, Taki FA, He Q, Liu H, et al. (2013) High-throughput deep sequencing shows that microRNAs play important roles in switchgrass responses to drought and salinity stress. Plant Biotechnology Journal 12:354–366.

31. Jinek M, Doudna JA (2009) A three-dimensional view of the molecular machinery of RNA interference. Nature 457:405–412.

32. Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB (2003) Prediction of mammalian microRNA targets. Cell 115:787–798.

33. Wu J, Zhu H, Song W, Li M, Liu C, et al. (2014) Identification of Conservative MicroRNAs in Saanen Dairy Goat Testis Through Deep Sequencing. Reproduction in Domestic Animals 49:32–40.

34. Sangokoya C, Telen MJ, Chi J-T (2010) microRNA miR-144 modulates oxidative stress tolerance and associates with anemia severity in sickle cell disease. Blood 116:4338–4348.

35. Kalimutho M, Blanco GDV, Di Cecilia S, Sileri P, Cretella M, et al. (2011) Differential expression of miR-144* as a novel fecal-based diagnostic marker for colorectal cancer. Journal of Gastroenterology 46:1391–1402.

36. Persengiev S, Kondova I, Otting N, Koeppen AH, Bontrop RE (2011) Genome-wide analysis of miRNA expression reveals a potential role for miR-144 in brain aging and spinocerebellar ataxia pathogenesis. Neurobiology of Aging 32:2316–e2317.

37. Rao PK, Missiaglia E, Shields L, Hyde G, Yuan B, et al. (2010) Distinct roles for miR-1 and miR-133a in the proliferation and differentiation of rhabdomyosarcoma cells. The FASEB Journal 24:3427–3437.

38. Chiba Y, Misawa M (2009) MicroRNAs and their therapeutic potential for human diseases: MiR-133a and bronchial smooth muscle hyperresponsiveness in asthma. Journal of Pharmacological Sciences 114:264–268.

39. He B, Xiao J, Ren A-J, Zhang Y-F, Zhang H, et al. (2011) Role of miR-1 and miR-133a in myocardial ischemic postconditioning. Journal of Biomedical Science 18:22.

40. Packer AN, Xing Y, Harper SQ, Jones L, Davidson BL (2008) The bifunctional microRNA miR-9/miR-9* regulates REST and CoREST and is downregulated in Huntington’s disease. The Journal of Neuroscience 28:14341–14346.

41. Laios A, O’Toole S, Flavin R, Martin C, Kelly L, et al. (2008) Potential role of miR-9 and miR-223 in recurrent ovarian cancer. Molecular Cancer 7:35.

42. Chen P, Price C, Li Z, Li Y, Cao D, et al. (2013) miR-9 is an essential oncogenic microRNA specifically overexpressed in mixed lineage leukemia-rearranged leukemia. Proceedings of the National Academy of Sciences 110:11511–11516.

43. Hudish LI, Blasky AJ, Appel B (2013) miR-219 Regulates Neural Precursor Differentiation by Direct Inhibition of Apical Par Polarity Proteins. Developmental Cell 27:387–398.

44. Yin D, Ogawa S, Kawamata N, Leiter A, Ham M, et al. (2013) miR-34a functions as a tumor suppressor modulating EGFR in glioblastoma multiforme. Oncogene 32:1155–1163.

45. Karaayvaz M, Zhai H, Ju J (2013) miR-129 promotes apoptosis and enhances chemosensitivity to 5-fluorouracil in colorectal cancer. Cell Death & Disease 4:e659.

46. Sakai A, Saitow F, Miyake N, Miyake K, Shimada T, et al. (2013) miR-7a alleviates the maintenance of neuropathic pain through regulation of neuronal excitability. Brain 136:2738–2750.