Identification and analysis of microRNAs in Botryococcus braunii using high-throughput sequencing

MicroRNAs (miRNAs) play important regulatory roles in the growth and development of organisms. The colonial green microalga Botryococcus braunii is an oil-rich algae and little is known about its miRNAs and their target genes. Here we constructed and sequenced a small RNA library for B. braunii using the HiSeq 2000 deep sequencing method. In total, we identified 42 known miRNA families and 14 novel miRNAs from B. braunii via sequence alignment and secondary structure prediction. Quantitative real-time PCR analysis indicated that the majority of miRNAs were effective and credible. Gene ontology analysis showed that most of the targets of miRNAs were classified as being involved in metabolic and cellular processes, gene expression regulation and stress/defense functions. Our findings provide the first large-scale identification and characterization of B. braunii miRNAs and their potential target genes. This study could lead to further identification of B. braunii miRNAs and enhance our understanding of their regulatory mechanisms in diverse biological and metabolic processes.

Microalgae have the potential to generate significant quantities of biomass and oil suitable for conversion to biodiesel.The colonial green microalga Botryococcus braunii has a remarkable ability to produce high levels of liquid hydrocarbons (Banerjee et al. 2002, Metzger & Largeau 2005).In previous studies, B. braunii was classified into 3 races (A, B and L) based on hydrocarbon structures (Kawachi et al. 2012, Yoshimura et al. 2013), and the genome and transcriptome of B. braunii were sequenced and analyzed (Weiss et al. 2011, Ioki et al. 2012a,b).To date, no systematic studies have been conducted on small RNA in this microalga.Thus, the aim of our study was to isolate miRNAs from B. braunii and (1) to establish a small RNA library of B. braunii, (2) to obtain sequences and information of miRNAs in B. braunii, and (3) to identify novel miRNAs in B. braunii using bioinformatics methods.The results will help us to identify the miRNA-based regulatory system of this microalga.

Strains and culture conditions
Botryococcus braunii was obtained from the Freshwater Algae Culture Collection, Institute of Hydro biology, Chinese Academy of Sciences, and propagated photoautotrophically in 1 l Erlenmeyer flasks containing 400 ml of modified Chu-13 medium (Largeau et al. 1980).The flasks were placed in a 25 ± 1°C illuminated incubator (Jiangnan Instrument Factory) for 14 d under a 14:10 h light:dark photoperiod and a light density of 50 ± 5 µE m −2 s −1 .The microalgal cells were then collected and cleaned with sterilized water.After being dried with hygroscopic filter paper, the samples were immediately frozen in liquid nitrogen and stored at −80°C before use.

Construction and sequencing of a small RNA library
Total RNA used for small RNA library construction was extracted from frozen B. braunii cells using Trizol (Invitrogen, Life Technologies) according to the manufacturer's instructions.Briefly, the frozen cells were re-suspended in 3 ml of Trizol reagent.Total RNA was then extracted and purified with the standard Trizol/ RNeasy column procedure.Following purification, total RNAs were digested with DNase and quantified by spectrophotometry.Total RNAs were then separated by size fractionation on denaturing 15% polyacrylamide gels.Fragments of 18 to 28 nt were gel-purified and then ligated to a 5'-adaptor and a 3'-adaptor.After reverse transcription, a 12cycle PCR reaction was performed, and the products were sequenced by the Beijing Genomics Institute (BGI) using the HiSeq 2000 deep sequencing method.

Small RNA analysis and miRNA identification
For analysis of small RNAs in B. braunii, unique small RNAs were compared and aligned with the sequences of non-coding RNAs (rRNA, tRNA, snRNA, snoRNA) available in Rfam (www.sanger.ac.uk/ software/Rfam) (Griffiths-Jones et al. 2005) and in the GenBank noncoding RNA database (www.ncbi.nlm.nih.gov/) using blastn with an e-value of 0.01 as the cutoff.In addition, all sequences were searched for conserved miRNAs of B. braunii using the plant mature miRNAs from miRBase (release 19.0) allow-ing for 2 mismatches and 3 gaps (Griffiths-Jones et al. 2008).

Prediction of novel miRNAs
Fold-back structures in miRNA precursors can be used to predict novel miRNAs (Allen et al. 2005).In this study, the Mireap program, developed by the BGI, was used to predict the novel miRNAs in B. braunii.The adopted strategy was as follows: (1) candidate miRNA sites were screened out from breakpoints defined by mapping of the small RNAs, (2) a minimal stringent criterion was used to select miRNA candidates and (3) the RNA secondary structure was checked with the Mfold program (Zuker 2003).

Prediction of miRNA targets
The miRanda program (www.microrna.org/)was used to detect potential target sites for the B. braunii candidate miRNA sequences.The parameters employed were as follows: match score S ≥ 90, target duplex free energy ΔG ≤ −20 kcal mol −1 and scaling parameter = 2.The miRNA-target duplexes were then checked manually according to rules suggested by John et al. (2004) and Huang et al. (2011).

Quantitative real-time PCR analysis
Ten miRNAs were randomly selected for quantitative real-time PCR (qRT-PCR) analysis.TaqMan miRNA assays were designed and ordered from Applied Biosystems.Each TaqMan miRNA assay includes a specific RT primer, forward and reverse primers and the TaqMan probe.Specific primers were designed according to the method of Chen et al. (2005) and synthesized by Sangon Biotechnology (Shanghai, China).Primers used in this study are listed in Table 1.qRT-PCR was performed according to the manufacturer's protocol and run on an ABI 7300 machine (Applied Biosystems) with thermal cycling parameters at 95°C for 10 min followed by 40 cycles of 95°C for 30 s, 60°C for 1 min according to the manufacturer's protocol.qRT-PCR reactions were run in triplicate with 2 biological replicates.The U6 snRNA was selected as a reference gene for normalization.A relative quantitative method (ΔΔC t ) was used to evaluate relative expression levels of different miRNAs in B. braunii.

Sequencing and analysis of the small RNA library
A small RNA library was established using Solexa technology to identify candidate miRNAs in Botryococcus braunii.In total, we obtained 10 020 143 sequence reads from the library.After filtering out low-quality tags, trimming adaptors and cleaning up shortages and contamination formed by adaptor− adaptor ligation, we obtained 9 786 842 (98.08%) clean reads, representing 2 005 344 unique sequences.Among the clean reads, 131 905 (1.35%) were similar to known miRNAs, and 6 602 059 (67.46%) were unannotated small RNA, suggesting that small RNA in B. braunii has not been investigated extensively in previous studies.Thus, our study has potential to discover more miRNA genes.The rest of the sequences were other types of RNA, including tRNA, rRNA, snRNA, snoRNA and other non-coding RNAs.The numbers and proportions of different categories of small RNAs are shown in Table 2.
As an important feature, size profile has often been used to distinguish miRNA from other small RNAs in previous studies.Most mature miRNAs with known functions commonly have a length of 20−24 nt (Wei et al. 2009).In this study, the length distribution pattern of the reads was analyzed and is presented in Fig. 1, showing that the majority of small RNAs in B. braunii library were 21 and 20 nt in size, accounting for 11.40 and 11.39% of the total reads, respectively (Fig. 1), followed by 19 nt (10.59%), 22 nt (10.17%), 23 nt (9.57%), and 24 nt (8.87%).This distribution pattern

Identification of known miRNAs in B. braunii
To identify the known miRNAs in B. braunii, the small RNA sequences were compared with the known plant miRNAs in miRBase 19.0 (www.mirbase.org/).Among the 9 786 842 sequences, 14 887 unique sequences were orthologs of known mi RNAs from other plant species previously deposi ted in miRBase (version 19.0).Allowing 1 or 2 mismatches between sequences, these miRNAs re presented 42 known miRNA families (Table 3).Solexa sequencing technology has the ability to generate millions of small RNA sequences, which could provide a resource with information on the abundance of various miRNA families and even distinguish between different members of a given family (Ruan et al. 2009).Thus, as a powerful tool, Solexa sequencing technology is often used to estimate the expression profiles of miRNA genes.In this study, the frequencies of known miRNA families in the sequenced library varied from 6 (miR164a) to 41 735 (miR408a-5p), indicating that expression varies significantly among different miRNA families.Among them, 13 miRNA families had more than 1000 sequence counts, suggesting these miRNAs were highly expressed in B. braunii.miR408a-5p, miR3948, and miR6426a were the most frequently expressed miRNA, with 41 735, 34 938, and 30 629 copies, respectively.Moreover, 15 miRNA families had 100 to 1000 sequence counts, and others had fewer than 100 sequence counts.2007).The predicted hairpin structures (Fig. 2) for the precursors of these novel miRNAs required 88 to 342 nt, and a majority of the identified miRNA precursors (71.4%) required 148 to 284 nt, which was more than what has been observed in other plants (Zhang et al. 2006).Moreover, the predicted novel miRNAs exhibited much lower expression levels.This finding was in line with the notion that the expression level of non-conserved miRNAs is lower than that of conserved miRNAs.Each novel miRNA family had only 1 member, and only the miR13 family had more than 1000 sequence reads.The majority of novel miRNA families had fewer than 20 reads.The low abundance of novel miRNAs might suggest a specific role for these mi -RNAs under various growth conditions, under biotic or abiotic stress, or during developmental stages (Zhao et al. 2010, Song et al. 2015).

Validation of novel miRNAs in B. braunii
To validate the predicted novel miRNAs, 10 mi -RNAs were randomly selected for qRT-PCR assays which were performed to examine whether the mi -RNAs were expressed in B. braunii.Fig. 3 shows that all 10 miRNAs were ex pressed in the samples of B. braunii, which was attributed to the conservation of miRNA families.The expression level of miR13 was relatively higher than that of the other 9 mi RNAs, which was in agreement with the Solexa sequencing data (Table 4).These findings were also in line with those reported by Gao et al. (2016), who found that miRNAs are relatively conserved in different species.Moreover, these results suggest that a majority of the miRNAs identified in B. braunii are effective and credible.In plants, miRNAs are involved in regulating a variety of biological processes, such as development, signal transduction, protein degradation, response to environmental stress and pathogen invasion (Lu et al. 2008).Their target sites have been shown to be primarily located in the coding regions.As expected, most miRNA target sites in B. braunii were located in the coding regions, and the putative target genes appear to be involved in a wide variety of biological processes (Fig. 4).Of these predicted targets, most were classified as being involved in metabolic and cellular processes, gene expression regulation, and stress/defense functions.Many metabolism networks have also been found, including lipid metabolism, amino acid metabolism, carbohydrate metabolism, energy metabolism, and nitrogen metabolism (Chi et al. 2011).

CONCLUSIONS
This is the first report of genome-wide identification of miRNAs in Botryococcus braunii using a highthroughput Solexa sequencing strategy.We have identified 42 known miRNA families and 14 novel miRNAs from B. braunii, suggesting that a significant number of novel microRNAs remain to be discovered and characterized.Additionally, we also predicted putative targets for these miRNAs.These results may help to improve our understanding of regulatory mechanisms of miRNAs in biological and metabolic processes of B. braunii.

Fig. 1 .
Fig.1.Lengths of unique small RNA sequences in Botryococcus braunii.The occurrences of each unique sequence read were counted to reflect the relative expression level, and only small RNA sequences in the range of 12 to 32 nt were considered

Fig. 3 .
Fig. 3. Quantitative real-time PCR (qRT-PCR) validation and expression analysis of miRNAs in Botryococcus braunii.Data are mean ± SD. qRT-PCR reactions were run in triplicate with 2 biological replicates

Table 3 .
The occurrences of each unique sequence read were counted to reflect the relative expression level, and only small RNA sequences in the range of 12 to 32 nt Putative miRNA families represented in the Botryococcus braunii miRNA pool