Identification and profiling of microRNAs of Euphausia superba using Illumina deep sequencing*

SONG Wei (宋煒), ZHU Yefei (朱葉飛), WANG Lumin (王魯民) , JIANG Keji (蔣科技) ,ZHANG Fengying (張鳳英), MA Chunyan (馬春艷), MA Lingbo (馬凌波) ,

1 Key Laboratory of Oceanic and Polar Fisheries, Ministry of Agriculture, East China Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Shanghai 200090, China

2 Key Laboratory of East China Sea Fishery Resources Exploitation, Ministry of Agriculture, East China Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Shanghai 200090, China

Abstract MicroRNAs (miRNAs) are an abundant class of conserved, non-coding small RNAs that play important role in gene regulation at post translational level. There have been no reports on the miRNAs of the Antarctic krill Euphausia superba despite the species’ crucial position in Antarctic food webs. Two small RNA libraries were constructed from eyestalk and muscle, subsequently, and deep sequencing analysis was performed to investigate and profile E. superba miRNAs. A total of 19 304 586 and 23 005 104 unique sequences were obtained from the eyestalk and muscle, respectively. After compared the small RNA sequences with the Rfam database, 12 342 039 and 7 907 477 reads in eyestalk and muscle were matched to the transcriptome sequence of E. superba. A total of 236 distinct miRNAs were identified after annotation to known animal miRNAs registered in miRBase 21. In both libraries, the most abundant known miRNA were miR-750 with 92 583 reads in muscle and miR-1304-3p with 56 386 reads in eyestalk while the average value was less than 106, revealing a wide range of different expression levels. In addition, miR-277a enriched in both libraries and may be involved in modulating Krebs cycle by targeting to Vimar. Differential expression analysis showed that 56 mature miRNAs were significantly up/down regulated according to expression fold change. It appeared that the expression of several abundant miRNAs maybe tissue-specific or tissue-bias. Notably, the expression pattern of miR-750 and miR-1 family, which was suggested as the crucial candidates, involved in muscle development. Taken together, this study provides the first miRNA profile of E. superba and some of these miRNAs were expected to play important role in immune response,reproduction, energy metabolism, and muscle development and so on. Thus, the results provides a reference for functional studies of miRNAs in E. superba.

Keyword: microRNA; Antarctic krill; Illumina deep-sequencing; eyestalk; muscle

1 INTRODUCTIO N

MicroRNAs (miRNAs) are an abundant class of conserved, approximately 18–25 nucleotides in length, endogenous non-coding RNA that regulate target gene expression through cleavage or translational inhibition (Cai et al., 2004), and then affecting crucial biological and metabolic processes.At least 35 828 miRNAs from hundreds of species,are recorded in the online repository miRBase 21(Kozomara and Griffiths-Jones, 2014). Though there is near-perfect pairing in plant miRNAs, animal miRNAs regulate target genes through imperfect sequence-specific binding to the 3′ or 5′ untranslated regions (Jopling et al., 2005) or even to the coding regions (Tay et al., 2008; Ryan et al., 2010), causing translational repression and, in certain cases,degradation of the mRNA (He and Hannon, 2004).Furthermore, regulation form of miRNA is not always inhibited, as the expression of tumor necrosis factor-alpha was up-regulated by miRNA-369-3p in human cells (Vasudevan et al., 2007). It is generally considered that miRNAs play the role of regulator at the post-transcriptional level. However, recent research has found that miR-24-1 activated enhancer RNA expression, altered histone modification, and increased the enrichment of p300 and RNA polymerase II at the enhancer locus (Xiao et al.,2017). Increasing evidence indicates that miRNAs are involved in embryogenesis and organogenesis(Melton et al., 2010; Pauli et al., 2011), cell development and migration, immune response(Brown et al., 2007; O’Connell et al., 2009), and hormone secretion (Cortez et al, 2011). Moreover,alternative expression of miRNAs can lead to series diseases or defective development (Calin et al, 2002).

High-throughput sequencing allows rapid and sensitive detection of miRNAs and can facilitate the discovery of novel, low abundance, and species- or tissue-specific miRNAs (Hafner et al., 2008). Using both deep sequencing and bioinformatics analysis to identify miRNAs in white shrimpLitopenaeus vannameiand giant freshwater prawnMacrobrachium rosenbergii, hundreds of conserved miRNAs were obtained (Tan et al., 2013; Xi et al., 2015). Antarctic krillEuphausiasuperbais possibly the world’s most abundant animal species, providing a massive protein resource that is targeted in the fisheries of numerous countries. However, with climate change and fishing operations, have resulted in a steady decline of the population density of Antarctic krill in the past decades (Atkinson et al., 2004; Cascella et al., 2015).For the sustainable development of this species,genomic and transcriptomic researches are expected to improve our understanding of the physiology and biology of krill.

This study aimed to identify and profile miRNA in the eyestalk and muscle tissues of Antarctic krill. The eyestalk is a complex and major organ for the secretion of hormones and involved in osmotic regulation,molting, epidermal color patterns, osmoregulation,modulation of glycaemia and reproduction (Turner et al., 2013; Katayama, 2016). The results will be useful to future research on the function of miRNAs involved in regulating development inE.superba.

2 MATERIAL AND METHOD

2.1 Sample collection

Healthy adult krillEuphausiasuperbafrom Antarctic waters of the Southern Ocean were collected on-board the Chinese icebreaking research vesselSnowDragonand then transported to the laboratory.The eyestalk and muscle tissues were promptly dissected, added to RNAfixer, and store at -80°C until RNA extraction.

2.2 Construction of small RNA library and Illumina deep sequencing

Total RNA was extracted using TRIzol?reagent(Invitrogen, CA, USA) following the manufacturer’s protocol. Total RNA quantity and purity were calculated with a 2100 Bioanalyzer and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number ＞7.0. Approximately 2 μg of purified total RNA was used to construct a small RNA library,according to the protocol of the NEBNext?Small RNA Library Prep Set for Illumina?(New England Biolabs, USA). Thereafter, the concentration and quality of small RNA librariy were measured using a Qubit?2.0 Fluorometer (Life Technologies) and an Agilent 2100 Bioanalyzer (Agilent Technologies).Finally, we performed the single-end sequencing(35 bp) on the Illumina HiSeq 2500 platform.

2.3 Bioinformatics analysis

Removing adaptors sequences, low quality tags,and sequences containing polyA tails allowed us to get clean reads and summarize the size distribution of small RNA. The retained 17–35 nt reads were mapped onto the reference transcriptome, filtering out unmapped sequences. The mapped reads were then classified into the tRNA, Repeat and Rfam database to remove rRNA, tRNA, snRNA snoRNA and other ncRNA, as well as the repeat sequence. Next, the sequences were mapped to intron and exon to identify degradation fragment of mRNA. After series modification, the retained sequences were blast to the miRNA precursor and mature miRNA of all animals in miRBase without mismatch, in order to identify conserved miRNAs.

MicroRNA read counts were normalized to the individual lane size by dividing each read count by the total number of reads in million per lane.Differentially expressed miRNAs between the two libraries were discovered using the R package DEGseq.P-values were calculated using DEGseq software. Scatter plots were generated based on log2 normalized read counts, with fold change = log2(normalized read counts in sample1/normalized read counts in sample2). Statistical significance was set at aP-value of ≤0.05 and fold change ≥2. To predict the targets of the miRNAs, we used sRNA-tools of target module.

Table 1 List of primers used for qPCR in this study

2.4 Quantitative Real Time-PCR

To validate and characterize the differentially expressed miRNAs identified in muscle and eyestalk of krill, we performed quantitative real-time PCR(qRT-PCR) busing the miRNA specific stem-loop primer approach (Table 1). Real-time PCR was performed using a standard SRYB Green PCR Kit protocol on an Applied Biosystems 7500 Sequence Detection System (P/N: 4329002, Applied Biosystems, Foster City, CA). The 20 μL PCR volume included 1 μL of RT product, 10 μL of 2X SYBR Green qPCR Master Mix, and 1 μL of primer(5 mmol/L each of the forward and reverse primers).The reactions were incubated at 95°C for 5 min,followed by 40 cycles of 95°C for 15 s, 60°C for 40 s.The relative expression of each miRNA was determined by comparativeCt(ΔΔCt) methods using U6 as the internal reference gene. Each experiment was repeated independently at least three times.

2.5 Statistical analysis

Statistical analysis was performed using SPSS v19.0 software. Data are presented as mean±standard error (S.E.). Statistical significance was determined using one-way analysis of variance (ANOVA) and post-hoc Duncan multiple range tests. Significance was set atP＜0.05.

Table 2 Summary of sequence counts from the two small RNA libraries

3 RESULT

3.1 General Features of the Solexa sequencing of small RNAs

This study generated two distinct small RNA libraries from the eyestalk (E) and muscle (M) to identify the Antarctic krill miRNAs using deep sequencing (Illumina Inc.). We obtained 21 609 615 and 24 755 663 raw reads from the E and M small RNA libraries, respectively. After filtering the junk reads, adapter sequences and the sequences outside the range of 18 to 35 nt, 19 304 586 reads from the E group and 23 005 104 reads from the M group were obtained for analysis (Table 2). Afterwards, we analyzed the length distribution of these clean reads to evaluate the sequencing quality. The length distribution of these clean sequence reads showed a peak at 20–24 nt for both libraries (Fig.1), consistent with the mainly products of Dicer derived cleavage(Yao et al., 2010).

Fig.1 Length distributions of the high quality filtered sequencing reads in the small RNA library of eyestalk and muscle for Antarctic krill

Fig.2 Annotation and classification of the small RNAs in the (a) eyestalk and (b) muscle library

Since the whole genomic sequence of the krill is unavailable, we analyzed the small RNAs from both libraries were analyzed using the de novo transcriptome from gill and hepatopancreas tissues ofE.superbaas the reference (unpublished data). Using Bowtie 2 v2.1.0 (Langmead and Salzberg 2012),63.93% of the clean reads for the E group and 34.37%for the M group could be mapped to the de novo transcriptome ofE.superbawithout any mismatches.Clean reads that could not be mapped to krill transcriptome, and matched sequence that did not form a hairpin structure, were removed. The mapped sequences were further categorized and annotated by aligning with sequences in Rfam and Repeat database(Fig.2). More than 80% of the mapping reads were identified as the unannotated (unano) group, while the smaller portion of the other RNAs lacked a genome resource. The repeat associated sequences from the E library were more various than those from the M library (Table 3).

3.2 Characterization of the miRNAs in eyestalk and muscle tissues of E. superba

To investigate conserved miRNAs in eyestalk (E)and muscle (M) ofE.superba, we aligned the mapped sequences with a non-redundant reference set of all animal miRNAs registered in miRBase 21 by allowing at most two mismatches outside of the seed region.The matched Small RNAs were identified as conserved miRNA orthologs inE.superba. We obtained 139 478 reads for E and 215 525 reads for M from the mapped reads as known miRNAs. Most of the predictedE.superbamature miRNA sequences,from 20 to 24 nt, for the E library started with an adenine residue at their 5′ end, which is consistent with findings for the aphidAcyrthosiphonpisum(Legeai et al., 2010) and the trematodeSchistosoma mansoni(De Souza Gomes et al., 2011), whereas the M library started with a thymidine residue. Ultimately,236 conserved mature miRNAs, which contained 64 miRNA-5p and 66 miRNA-3p, were identified from the two libraries, and 48 of them were co-expressed in both libraries.

Table 3 Repeat sequences from the small RNA libraries

After normalizing the individual mapped reads, the most abundant miRNAs identified in each library were listed in in descending order of abundance(Tables 4, 5); these results indicated that the expression of several of the abundant miRNAs may be tissue- specific. The most abundant of the known miRNA was miR-750, with 92 583 reads detected in muscle,with an average value of 68. Among the small RNA library of eyestalk, miR-1304-3p was the most accumulated of the known miRNAs, and showed 56 386 reads, with an average value of 105. Especially,miR-277a was abundantly expressed, with more than 10 000 reads in both libraries.

Table 4 Most abundant mature miRNA detected in the eyestalk library

Table 5 Most abundant mature miRNA detected in the muscle library

Fig.3 Differential expression of known miRNAs between eyestalk and muscle

3.3 Differential expression of miRNAs between eyestalk and muscle

We compared the miRNA expression levels between the two small RNA libraries of eyestalk and muscle tissues. For this purpose, 56 mature miRNAs were significantly up/down regulated in terms of expression fold change, based on a thresholds of 2-times fold change (FC) and an overallP-value of＜0.05. Among the most abundant miRNAs that were expressed in relatively high amounts in the two libraries (Fig.3), miR-1304-3p (-2 log2 FC), miR-2478 and miR-6561-5p (both -1.5 log2 FC), were significantly down regulated, whereas miR-750, miR-1b and miR-8-3p (all ～7 log2 FC), were significantly up regulated compared to eyestalk. Interestingly, the expression of three miR-1 members (miR-1a, -1b,and -1-3p) was significantly higher in muscle than in eyestalk, showing muscle-biased expression.

3.4 Quantitative real-time PCR validation

To confirm the accuracy of the sequencing results,four miRNAs were selected and determined by qRTPCR via miRNA specific stem-loop primers (Fig.4).Among the four miRNAs, miR-1a and miR-1304-3p were relatively high expressed in muscle and eyestalk respectively, and miR-317-3p and miR-184 were coexpressed in the two organs. The relative expression levels of these miRNAs were consistent with the results of Solexa sequencing.

4 DISCUSSION

The increasing evidence demonstrates that miRNAs as an abundant group of fundamental regulators, participating in many biological and phsiological processes (Begemann, 2008). In the last decades, over thirty thousand of miRNA have been identified in various species (miRBase v21).Improvements in RNA deep sequencing technologies and bioinformatics have improved detection of miRNA transcriptome detection in crustacean, which are lacking a reference genome, such as in freshwater crayfishProcambarusclarkia(Ou et al., 2013), mud crabScyllaparamamosain(Li et al., 2013), mitten crabEriocheirsinensis(Song et al., 2014), and tadpole shrimpTriopscancriformis(Ikeda et al.,2015). However, no data are available regarding miRNAs of krill. In this study, were constructed two distinct small RNA libraries of from the eyestalk (E)and muscle (M) tissues, to identify the Antarctic krill miRNAs by high-throughput sequencing technologies.We obtained 19 304 586 and 23 005 104 clean reads from the small RNA libraries of E and M, respectively.The length distribution analysis revealed a peak at the range of 20–24 nt for both libraries (Fig.1), which matched to the typical length products for Dicer derived cleavage (Yao et al., 2010). Using the mRNA,Rfam and Repeat databases, more than 80% of the mapped reads were classified as the unannotated(unano) group, and the other small RNAs accounted for the smaller proportion. These results suggested the limitations of the reference transcriptome and the high quality of these new libraries.

Fig.4 Comparison of the expression levels of four miRNAs in eyestalk (E), muscle (M) and hepatopancreas (H) discovered through qRT-PCR

In this study, a total of 236 conserved miRNAs from eyestalk and muscle tissues were identified and annotated. We then focused on the most abundantly expressed and conserved miRNAs in the theE.superbalibraries (Tables 3, 4). Additional sequencing is required for identification of low expression level and novel miRNAs. The most abundant known miRNAs in the libraries were miR-750, with 92 583 reads in muscle, and miR-1304-3p,with 56 386 reads in eyestalk, while the average values were less than 106, revealing a wide range of different expression levels. The polymorphic distribution of miRNAs often suggests a role in distinct tissues or biological processes. As shown in previous research, miR-317, miR-2, miR-1175 and miR-277 are evolutionary conserved miRNAs, having similar sequences and functions among species (Ruby et al., 2007). The miR-317 is involved in sex peptides reverse (Fricke et al., 2014) and and possibly targetsDrosophilacyclin B during the G2-M transition(Pushpavalli et al., 2014). The miR-2 demonstrated high expression in the first metaphase (MI) of meiosis and luciferase reporter gene assay analysis showed that miR-2 can down regulate the crab cyclin B gene(Song et al., 2014). In addition, miR-2 has a potential function in neurodevelopment in theDrosophila(Marco et al., 2012). The miR-1175 is regarded as key regulator of quinone oxidoreductase and glutathione peroxidase, which are important to the immune system (Ou et al., 2013). The expression level of miR-277-3p in the goldenrod gall flyEurostasolidaginiswas significantly reduced in frozen larvae (Courteau et al., 2012). The Vimar enzyme is employed as an indicator of mitochondrial oxidative capacity and was revealed to be under the regulation of miR-277(Bruce, et al., 2004; Bos, 2005). In this study, maybe due to the environment change during the krill catches, miR-277a was found to be a dominant miRNA, detected in both muscle and eyestalk, which suggests that miR-277 might be involved in freeze tolerance and energy metabolism in the Antarctic krill.

We further compared the transcriptional levels of the miRNAs between the two small RNA libraries of eyestalk and muscle tissues. In terms of differently expressed miRNAs between the libraries, 34 miRNAs were up regulated and 22 miRNAs were down regulated in muscle as compared with eyestalk.Notably, miR-750 was the most abundant miRNA of this study, by more than 7-fold in muscle as compared to eyestalk. In addition, all the three members of the miR-1 family showed muscle-biased expression. As demonstrated by previous study (Huang et al., 2014),miRNAs mediate cell proliferation and differentiation of muscles by modulating muscle development related gene expression. Conditional knockout (cKO)of Dicer1 (a distinct protein that participate in the biosynthesis of miRNAs) from skeletal muscle leaded to myofibers hypoplasia and apoptosis in mouse(O’Rourke et al., 2007). It was found that miR-1 is restrictive or abundantly expressed in skeletal muscle tissue, which suggests an important role in regulating muscle development and growth (Yan et al., 2012).Moreover, in zebrafish, actin organization and sarcomere assembly were disrupted under defective miR-1 expression (Mishima et al., 2009). Analogously,loss-function of miR-1 inDrosophilalead to serious amyotrophy (Kwon et al., 2005; Sokol and Ambros,2005). MHCs (myosin heavy chains), a major component of myosins, were identified as a potential growth marker gene (Harrington and Rodgers, 1984).In the swimming crabPortunustrituberculatus, miR-750 has a higher expression level in larger-sized individuals than in smaller individuals, and MHC was the only predicted target gene of miR-750, indicating that it may stimulate growth by regulating MHC expression (Ren et al., 2016). Overall, the miR-750 and miR-1family are the key candidates involved in krill muscle development.

We found that the two most abundant miRNAs(miR-1304-3p and miR-6561-5p) were significantly up regulated in eyestalk as compared with those in muscle. Intriguingly, both of those miRNAs were also new candidates detected in embryonic stem cells(Morin et al., 2008; Shao et al., 2012), signifying their involvement in cell differentiation, migration and apoptosis., although the available information on miRNAs in crustaceans is quickly expanding, most research have focused on the immune response (Li et al., 2013), responses to viral infection (Li et al., 2013)and stress responses (Lv et al., 2016), as detected mainly in hepatopancreas, gills and gonad tissues.However, few studies have considered the expression patterns and functions of miRNAs in eyestalk of crustacean. This is the first profiling of miRNAs in krill, and the research has identified hundreds of miRNAs in eyestalk and muscle, some of these miRNAs displaing a potential role in the immune response, reproduction, energy metabolism and muscle development. These data provide a theoretical basis for the functional research of miRNA inE.superba.

5 CONCLUSION

To our knowledge, this is the first study on miRNA expression profiles of eyestalk and muscle tissues in krill. A total of 236 distinct miRNAs were identified in the small RNA libraries from eyestalk and muscle of the krill. Differential expression analysis demonstrates that some of the identified miRNAs are tissue-specific or tissue-biased, which suggested potential roles in the regulation of endocrine function,growth and development of krill. Specially, miR-750 and miR-1 were found to be highly expressed in krill muscle, revealing a significant role in krill muscle development. The results presented here will be of assistance in miRNAs functional studies ofE.superba, such as aspects of the species’ in immune response, reproduction, energy metabolism, and muscle development.

6 DATA AVAILABILITY STATEMENT

The datasets during and/or analyzed during the current study are available from the corresponding author on reasonable request.

7 ACKNOWLEDGMENT

We thank Liwen Bianji, Edanz Group China (www.liwenbianji.cn), for editing the English text of a draft of this manuscript.