High throughput sequencing of small RNAs transcriptomes in two Crassostrea oysters identifies microRNAs involved in osmotic stress response.

Increasing evidence suggests that microRNAs post-transcriptionally regulate gene expression and are involved in responses to biotic and abiotic stress. However, the role of miRNAs involved in osmotic plasticity remains largely unknown in marine bivalves. In the present study, we performed low salinity challenge with two Crassostrea species (C. gigas and C. hongkongensis), and conducted high-throughput sequencing of four small RNA libraries constructed from the gill tissues. A total of 202 and 87 miRNAs were identified from C. gigas and C. hongkongensis, respectively. Six miRNAs in C. gigas and two in C. hongkongensis were differentially expressed in response to osmotic stress. The expression profiles of these eight miRNAs were validated by qRT-PCR. Based on GO enrichment and KEGG pathway analysis, genes associated with microtubule-based process and cellular component movement were enriched in both species. In addition, five miRNA-mRNA interaction pairs that showed opposite expression patterns were identified in the C. hongkongensis, Differential expression analysis identified the miRNAs that play important regulatory roles in response to low salinity stress, providing insights into molecular mechanisms that are essential for salinity tolerance in marine bivalves.

MicroRNAs (miRNAs) are a class of endogenous small noncoding RNAs, which regulate gene expression post-transcriptionally in animals and plants1. MiRNAs are abundant in animal genomes and have been reported to play primary roles in a broad range of biological and metabolic processes, such as regulation of individual development23, abiotic stress response45, pathogen defense and innate immune response678. In a previous study, more than 60% of human protein-coding genes are predicted as targets of miRNAs based on in silico analysis9.

Recent progress in the development of genomic techniques, especially high-throughput sequencing, has greatly facilitated transcriptome analysis of ecologically and economically important animals. Next-generation sequencing based RNA-seq analysis has been widely used to uncover expression patterns under different conditions1011. Recently, high-throughput sequencing of small RNA transcriptomes for miRNA discovery and expression profiling has been conducted in numerous model organisms1213 and plants414. Next-generation sequencing technologies have also been used for miRNA studies in a number of non-model organisms including several fish species1516171819, echinoderms2021222324, and marine bivalves72526272829.

Crassostrea is one of the most studied bivalve genera due to its worldwide distribution, strong adaptability and high economic importance. The ability to cope with abiotic and biotic stresses is vital to survival of the oysters because of their intertidal inhabiting lifestyle. It has long been an interesting question how the oysters could be survived in response to the environmental stresses involved in the fluctuation of salinity, and air exposure. It is reported that the biogenesis of miRNA, the expression of mRNA targets, and the activities of miRNA-protein complexes were actively regulated under stress conditions30. Many studies have been performed to reveal the alterations of mRNA expression313233343536 and protein expression3738394041 under abiotic and biotic stresses in the oyster42, while expression profiling of miRNAs under osmotic stress remains largely unexplored.

The Pacific oyster (C. gigas) and the Hongkong oyster (C. hongkongensis) are two major Crassostrea species along the coasts of China43. These two oyster species form distinct fauna assemblages. C. gigas can live in the environments with salinity over 20‰44, while C. hongkongensis inhabits in the estuarine condition with river disemboguing4546. Due to the drastic differences in salinity tolerance, the C. gigas and C. hongkongensis provide an excellent model system to study the roles of miRNAs in adaptation to hypo-osmotic stress.

In this work, we performed high-throughput sequencing of small RNA transcriptomes of gill tissues in C. gigas and C. hongkongensis after osmotic stress treatment. The gill tissue is a primary interface between the hemolymph or cytoplasm and the external environment in marine molluscs whose osmolality fluctuates widely47. Gill is known to be responsive to environmental challenges at the transcript level3548 and proteomic level41. Based on the recently released genome sequence of the Pacific oyster49 and the global expression profile of mRNA transcripts in the gill tissues from our previous study48, we identified known and novel miRNAs from the two oysters. The integrated analysis of miRNA and mRNA expression profiles in the two oyster species allowed identification of miRNA-mRNA interaction pairs. Further Gene Ontology and pathway analysis enabled the investigation of their putative biological functions. This work generated genomic resources of miRNAs that are valuable for further studies in the oysters, and provided insights into the roles of miRNAs in response to low salinity stress.

Results

High-throughput sequencing of small RNAs

Four small RNA libraries constructed from the two groups of C. gigas (one in salinity 25‰, PC; one in salinity 8‰, PT) and two groups of C. hongkongensis (one in salinity 25‰, HC; one in salinity 8‰, HT). High-throughput small RNA transcriptome sequencing generated a total of 13,268,157, 14,682,519, 22,288,517, and 21,645,549 raw reads from PC, PT, HC and HT groups, respectively. After removing low quality reads and adaptors, a total of 320,436, 429,972, 464,695 and 572,293 unique clean reads from PC, PT, HC and HT groups, respectively. These unique clean reads were further annotated based on the Pacific oyster genome, and were classified into miRNAs, tRNAs, rRNAs and others by blasting with Rfam and miRBase (Table 1). The reads annotated with other noncoding RNAs rather than miRNAs were excluded for further analysis. The unannotated RNAs were used to identify novel miRNAs.

Table 1

Summary of small RNA transcriptome sequencing of gill tissues in C. gigas and C. hongkongensis.

Category	C. gigas				C. hongkongensis
	Control		Treated		Control		Treated
	Unique	Redundant	Unique	Redundant	Unique	Redundant	Unique	Redundant
miRNA	1,878	9,150,470 (71.38%)	2,108	9,492,885 (67.12%)	1,223	14,682,335 (68.04%)	1,401	13,579,598 (65.01%)
rRNA*	5,968	86,221 (0.67%)	8,297	146,244 (1.03%)	7,552	143,791 (0.67%)	9,732	213,415 (1.02%)
tRNA*	976	31,809 (0.25%)	1,491	88,173 (0.62%)	136	1,978 (0.009%)	220	3,367 (002%)
Other Rfam category*	1,888	6,843 (0.05%)	2,472	12,042 (0.09%)	2,396	9,891 (0.04%)	2,832	15,273 (0.07%)
Unannotated	56,073	823,478 (6.42%)	68,386	1,081,813 (7.65%)	103,777	892,043 (4.13%)	121,226	1,169,750 (5.60%)
Mapping to genome*	66,783	10,098,821 (78.78%)	82,754	10,821,155 (76.51%)	115,084	15,730,038 (72.89%)	135,411	14,981,403 (71.72%)
Clean reads	320,436	12,819,335	429,972	14,142,838	464,695	21,580,148	572,293	20,889,308
Raw data	13,268,157		14,682,519		22,288,517		21,645,549

*excluded for further analysis.

The size distribution of all miRNAs indicated that the majority of miRNAs were with lengths of 21–23 nt long (Fig. 1). The majority of reads were with length of 22 nucleotides in the four libraries miRNAs occupied the major percentages (52.7% for PC, 56.1% for PT, 52.9% for HC, 56.4% for HT) in the four libraries.

Figure 1

Length distribution of small RNAs in four groups from C. gigas and C. hongkongensis.

Identification of conserved miRNAs

A total of 137 and 85 conserved miRNAs were identified in the C. gigas (PC: 132, PT: 137) and C. hongkongensis (HC: 82, HT: 81), respectively. One hundred and seven miRNAs in C. gigas and 54 miRNAs in C. hongkongensis with read counts greater than 10 were listed in Table S1. A total of 39 miRNAs were identified in both species.

Comparison of the expression of all conserved miRNAs within the two oyster species revealed that miR-10a was expressed with the highest level in each group of both species. The expression of miR-10a was followed by miR-184, miR-184-3p, miR-10b and miR-981 in C. gigas, while was followed by miR-981, miR-8, miR-10-5p and miR-2001 in C. hongkongensis (Table S1).

Identification of novel miRNAs

The precursors of several unannotated miRNAs were identified based on the C. gigas genome49 and C. hongkongensis transcriptome48. A total of 65 novel miRNAs were identified in the C. gigas and two novel miRNAs were identified in the C. hongkongensis, respectively (Table S2). Of the 65 novel miRNAs from C. gigas, 32 and 37 novel miRNAs were identified in PC and PT groups, respectively, with four novel miRNAs being identified in both groups. In C. hongkongensis, only one novel miRNA was identified in both groups and the other one was only identified in the HT group.

Identification of differentially expressed miRNAs after osmotic stress

Based on the statistical analysis of read counts in the four miRNA libraries, a total of six differentially expressed miRNAs (three up-regulated and three down-regulated miRNAs) between PC and PT groups were identified in the C. gigas and two differentially expressed miRNAs (one up-regulated and one down-regulated miRNAs) between HC and HT groups were identified in the C. hongkongensis (Table 2). The miR-2353 was identified to be down-regulated in both species.

Table 2

Differentially expressed miRNAs identified from the C. gigas and C. hongkongensis between control and treatment groups.

miRNA	Read counts (Control group)	Normalized read counts (Control group)	Reads_counts (Treatment group)	Normalized read counts (Treatment group)	Log2(foldchange)	P-value	Species
scaffold43364_10952	23	23.62598138	588	572.420666	4.598629596	3.74E-05	C. gigas
cgi-miR-2353	72	73.6333131	2	1.955636572	−5.234648421	0.00174025	C. gigas
cgi-miR-1984	164,760	169,244.204	534,981	520,806.4292	1.62164083	0.00184559	C. gigas
cgi-miR-183	2,160	2,208.999393	908	887.8590035	−1.31499053	0.02035668	C. gigas
cgi-miR-92-3p	73,928	75,605.0496	165,791	162,113.4714	1.100449483	0.02711921	C. gigas
cgi-miR-184-3p	601,642	618,016.6388	308,210	300,043.8325	−1.042472405	0.04226714	C. gigas
chk-miR-3205	2	1.992031841	110	110.4400018	5.792878298	8.05E-09	C. hongkongensis
chk-miR-2353	118	117.5298786	4	4.016000064	−4.871124465	3.85E-08	C. hongkongensis

Validation of differentiallly expressed miRNAs by qRT-PCR

Those six differentially expressed miRNAs (scaffold43364_10952, cgi-miR-1984, cgi-miR-92, cgi-miR-183, cgi-miR-2353, cgi-miR-184-3p) identified in C. gigas and two differentially expressed miRNAs (chk-miR-3205 and chk-miR-2353) identified in the C. hongkongensis were validated by quantitative real-time PCR. All of these miRNAs showed a consistent expression pattern with the results from small RNA sequencing (Fig. 2), indicating high reliability of the analysis.

Figure 2

Expression of eight miRNAs determined by qRT-PCR.

The eight miRNAs included scaffold43364_10952 (a), miR-92-3p (b), miR-1984 (c), miR-183 (d), miR-184-3p (e), and miR-2353 (f) in C. gigas and miR-2353 (g) and miR-3205 (h) in C. hongkongensis. 5S gene was used as an internal control to calibrate the cDNA template for all the samples. Each values were shown as mean ± SD (n = 6).

Prediction and annotation of miRNA target genes

To better understand the functions of the identified miRNAs, candidate target genes were predicted using the genome of C. gigas and the transcriptome of C. hongkongensis as the reference transcript set by the miRanda software. In C. gigas, 472 target genes were predicted for the three up-regulated miRNAs and 880 target genes were identified for the three down-regulated miRNAs. In C. hongkongensis, 739 target genes were identified for the up-regulated miRNA and 3,607 target genes were identified for the down-regulated miRNAs. The GO assignment distribution of target genes was shown in Fig. 3. The biological functions of these target genes were further investigated using KEGG pathway analysis. More than 200 different pathways were found, and the most frequently represented pathways were involved in signal transduction and immune system (Fig. S1). This annotation of genes targeted by differentially expressed miRNAs suggests that the gene functions regulated by these miRNAs in response to how salinity stress are similar in the two species. The results indicate that there is no significant difference in gene function regulated by miRNAs in response to low salinity stress between the two oyster species.

Figure 3

GO distribution of target genes of differentially expressed miRNAs of C. gigas (a) and C. hongkongensis (b). Red: target genes of up-regulated miRNAs; Purple: target genes of down-regulated miRNAs.

The basic biological function of each putative target gene was classified based on the Gene Ontology enrichment and KEGG pathways analysis. The hypergeometric test was used to identify the overrepresented GO terms and KEGG pathways.

The GO enrichment analysis provided several important biological processes enriched in both C. gigas and C. hongkongensis, such as microtubule-based process, cellular component movement, catabolic and metabolic process of purine nucleoside, and intracellular signal transduction. The enriched GO terms shared between the two species were shown in Fig. 4 and listed in Table S3. All the enriched GO terms were listed in Table S4 (for C. gigas) and Table S5 (for C. hongkongensis). Furthermore, the pathway analysis suggested that ECM-receptor interaction played an important regulatory role in both species. Several pathways were specifically enriched in the two species, such that cell adhesion molecules was enriched in the C. gigas, while MAPK signaling pathway and Focal adhesion were enriched in the C. hongkongensis (Table 3).

Figure 4

Venn diagram of shared enriched GO terms between C. gigas and C. hongkongensis.

Table 3

List of significantly enriched KEGG pathways of candidate targets of differentially expressed miRNAs in C. gigas and C. hongkongensis.

KEGG pathway	Term	Gene Counts	Gene number of Pathway	P-value	species
Map04514	Cell adhesion molecules (CAMs)	6	12	0.000102	C. gigas
Map04512	ECM-receptor interaction	7	27	0.002859	C. gigas
Map04010	MAPK signaling pathway	30	113	0.000105	C. hongkongensis
Map04510	Focal adhesion	24	85	0.000185	C. hongkongensis
Map04515	Rap1 signaling pathway	20	79	0.002729	C. hongkongensis
Map04512	ECM-receptor interaction	8	21	0.003862	C. hongkongensis
Map00410	Beta-Alanine metabolism	9	27	0.006183	C. hongkongensis
Map05164	Influenza A	17	69	0.007415	C. hongkongensis
Map04062	Chemokine signaling pathway	14	54	0.009128	C. hongkongensis

Integrated expression analysis of microRNAs and their target mRNAs during osmotic stress

Gene expression regulation by miRNA relies on the miRNA-mRNA interactions. Based on the fact that miRNA regulate gene expression by inhibiting translation or inducing deadenylation of mRNA followed by their degradation5051, expression profiling of miRNAs and mRNAs should reveal an inverse relationship if one regulates the other. Here we compared the predicted target genes of differentially expressed miRNAs in the two Crassostrea species by miRanda with differentially expressed mRNAs of the same samples48. As a result, expressions of two miRNAs, including chk-miR-2353 and chk-miR-3205, were negatively correlated with that of their corresponding target genes in response to low salinity stress (Table 4).

Table 4

Interaction pairs of differentially expressed miRNAs and predicted mRNAs in C. hongkongensis.

miRNA name	Target Gene-Interaction	miRNA expression level		mRNA expression level		Target annotation
		log2(fold change)	P-value	log2(fold change)	P-value
chk-miRNA-2353	comp106174	−4.87	3.85E-08	5.75	9.30E-06	ATP grasp domain containing protein 1
	comp24781	−4.87	3.85E-08	6.05	3.15E-05	cAMP responsive element binding protein like 2
chk-miRNA-3205	comp41123	5.79	8.05E-09	HC	6.90E-04	Hygromycin phosphotransferase
	comp54488	5.79	8.05E-09	HC	2.38E-05	Replication factor A protein 1
	comp62265	5.79	8.05E-09	−6.40	1.44E-05	von Willebrand factor D and EGF domain containing protein

HC: mRNA expressed only in HC group with no value of log2 (fold change).

Of the osmoregulatory candidate genes identified in C. gigas49, 103 were identified in the target genes of differentially expressed miRNAs of C. gigas (Table S6) and 570 in those of C. hongkongensis based on 1237 transcripts of C. hongkongensis blasted to the osmoregulatory candidate genes of C. gigas (Table S7).

Discussion

Since their roles in post-transcriptional regulation were unraveled, miRNAs have been extensively investigated in many organisms. Although little progress has been made in non-model species, the advances of high throughput sequencing technology provide unprecedented opportunities to efficiently characterize small RNA transcriptome in the molluscs. In this study, we conducted a high-throughput small RNA transcriptome sequencing in two oyster species in order to determine the roles that miRNA play in response to low salinity stress. For the consistent treatment to the two Crassostrea species, salinity 25‰ was considered as an optimal salinity level. The salinity level at their original habitat (the collection site) was 25‰ fluctuating with tide frequently. In previous studies, salinity 25‰ was used as the salinity of acclimation and proved as an appropriate salinity for physiological activities in both oyster species45525354. Salinity 8‰ is almost a extreme salinity for C. gigas4155, which was used as treated salinity level for osmotic stress.

The small RNA transcriptome sequencing resulted in identification of 202 and 87 miRNAs from C. gigas and C. hongkongensis, respectively. There is a significant difference in the number of miRNAs identified in the two oyster species. It may be mainly due to the significant difference in abundance of the genomic information of the two oyster species. Further expression analysis indicated that six miRNAs in C. gigas and two miRNAs in C. hongkongensis were potentially involved in regulating the acclimation to low salinity stress conditions.

High-throughput sequencing provided sufficient and reliable genomic data for downstream small RNA analysis. The length distribution of reads (Fig. 1) in the four libraries was consistent with the observation in other aquatic organisms202556, suggesting the conservation of miRNAs and the high quality data obtained for downstream analysis in this work. Using the genome of the Pacific oyster and the transcriptome of the Hongkong oyster as references for miRNA mapping, 202 miRNAs were identified in the gill of the C. gigas, with 137 being known miRNAs and 65 being novel miRNAs. Of the 137 known miRNAs, 71 known miRNAs were reported in a previous study that were identified from the haemocytes of C. gigas25. The remaining 66 known miRNAs could be a reflection of different miRNAs expression between the gill and haemocyte of 65 novel miRNAs, 16 miRNAs were also reported in haemocytes of the Pacific oyster25, and the remaining 49 novel miRNAs were discovered only in this study.

The expression of miRNAs under normal conditions was consistent with observations in other organisms. The miR-10 family (miR-10a, miR-10b, miR-10-5p), miR-184 family (miR-184, miR-184-3p), and miR-981 were the most abundant ones in the control groups of the two oyster species. Enhanced expression of miR-10a and miR-184 was also observed higher in haemocytes of C. gigas25, sea cucumber (Apostichopus japonicas)2021 and Chinese mitten crab57. These miRNAs were conserved and had a long evolution history from Nephrozoa58, suggesting that they might play crucial roles in essential biological processes in various organisms596061.

Cgi-miR-183, cgi-miR-184-3p, cgi-miR-2353, cgi-miR-1984, cgi-miR-92 and scaffold43364_10952 were differentially expressed in the gills of C. gigas under osmotic stress, while chk-miR-3205 and chk-miR-2353 were differentially expressed in the gills of C. hongkongensis. Previous studies reported that miR-183 was involved in response to environmental stress. The miR-183 functioned as a regulator of target genes to activate the c-Jun N-terminal kinase mitogen-activated protein kinase pathway and regulate the pathways of apoptosis6263. Moreover, members of miR-183 family were found to be expressed in innervated regions of invertebrate deuterostomes64 and regulate down-stream effectors functioning in actin cytoskeleton and plasma membrane65. It is reported that cytoskeleton rearrangement may maintain internal and external osmotic pressure balance when the Pacific oyster responds quickly to osmotic stress34. In this study, cgi-miR-183 was significantly down-regulated under low osmotic conditions compared to normal conditions, which may regulate the genes involved in cytoskeleton to response to salinity stress. The miR-184-3p is a specific miRNA that can modify the microRNP (microribonucleoproteins) function and relieve the repression induced by stress66. The miR-92, belonging to miR-17-19 cluster, is an immune-related miRNA67, and is associated with pro-proliferative and anti-apoptotic properties68. Notably, the differentially expressed miRNAs with known functions identified in the present study are mainly involved in immune-related function and response to abiotic stress. However, the function of cgi-miR-1984, miR-2353, chk-miR-3205 and scaffold43364_10952 are still unknown. The miR-1984 was first discovered as a gastropod-specific miRNA gene58 and was only found in Lottia gigantea and Haliotis rufescens. However, the identification of miR-1984 in the haemocytes of C. gigas in the previous study25 and gills of C. gigas in this study suggested that miR-1984 was not gastropod-specific but could be mollusk-specific. The expression of miR-1984 increased significantly in response to low salinity stress, heat stress and bacteria challenge suggesting that it could be implicated in certain physiological functions such as oxidation reduction and energy metabolism. The miR-2353 and miR-3205 were previously discovered in cattle69 and silkworm70. The target genes of the three miRNAs (miR-1984, miR-2353, and miR-3205) might be implicated in signal transduction and proteolysis for energy metabolism and other physiological functions. Further functional analysis are needed to elucidate their roles in molluscs.

Integrated expression analysis of miRNAs and their target genes identified the interaction of miRNAs and mRNAs involved in regulating specific biological processes. The enrichment analysis of putative target genes revealed the relationship between miRNAs expression and response to low-osmotic stress in the two oyster species. Based on the results of differentially expressed transcripts in C. gigas and C. hongkongensis48, we were able to identity a total of five functional miRNA-mRNA interaction pairs in response to low salinity stress (Table 4). As anticipated, all identified miRNA-mRNA interactions showed that the expression of miRNAs had negatively correlationships with that of their mRNA targets, consistent with the observations that miRNAs regulate target gene expression by repressing their targets71 through transcript cleavage or translation repression.

The chk-miRNA-2353 appeared to target ATP grasp domain containing protein 1 (ATPGD1) and cAMP responsive element binding protein like 2 (CREBL2). ATPGD1 is involved in metabolic process and catalyzes the degradation of beta-alanine in KEGG to maintain osmotic equilibrium under hypo-osmotic stress in the oysters. Its expression was found to increase significantly in the C. gigas on the 7th day after hypo-osmotic stress, and reached the highest level under the condition with salinity of 1055. The CREB family was reported to be involved in immune processes in mollusca7273. The chk-miRNA-3205 targets hygromycin phosphotransferase, Replication factor A protein 1 (RFP1) and von Willebrand factor D and EGF domain containing protein (VWDE). RFP1 has zinc finger domain, which has been reported to be involved in stress response and apoptosis in previous studies347475. VWDE was putatively related to immunity and stress7677. These miRNA-mRNA interaction pairs deserve future investigation.

Dramatic differences in expression profiles of miRNAs were observed between the two oyster species. We detected six miRNAs and 48 corresponding target genes in the C. gigas, while two miRNAs and 408 corresponding target genes were identified in the C. hongkongensis, in response to hypo-osmotic stress. Five miRNA-mRNA interaction pairs were found in C. hongkongensis, while none was identified in C. gigas. This result was far below observation reported in previous studies7578. Possible reasons to explain this observation include: (1) the computationally-predicted targets do not represent the actually existing interactions, (2) strict filtering criterion may exclude certain actual interaction pairs, (3) the expression of specific genes are regulated at temporal and spatial levels, and most targets may not be regulated at this point, and (4) there is less interaction pairs in the C. gigas than that in the C. hongkongensis, which is consistent with the observation that C. hongkongensis had higher levels of tolerance than that of C. gigas in response to acute hypo-osmotic stress. Similarly, comparison between the control group of C. gigas and C. hongkongensis also showed a dramatic difference, indicating the different tolerance to low salinity stress between the two oyster species. Moreover, miR-2353 may be the main miRNAs to resist low osmotic stress, which was the only one miRNA down-regulated in both species of the osmoregulatory candidate genes identified in C. gigas49, the targets of differentially expressed miRNAs in C. gigas and C. hongkongensis that were identified in osmoregulatory candidates may be the candidate genes regulated by miRNAs in response to osmotic stress. There were more osmoregulatory candidates regulated by miRNA in C. hongkongensis than in C. gigas. The results indicated that C. hongkongensis can regulate more osmoregulatory genes in response to acute hypo-osmotic stress.

Materials and Methods

Ethics statement

The C. hongkongensis and C. gigas are not endangered or protected species, and there is no requirement for permission to perform experiments involving the two oysters (invertebrate) in China.

Sample collection

The samples used in this work were same as described in our previous RNA-seq study48. In brief, adult individuals of C. gigas (2 year-old) were collected from Weihai (Shandong, China), and adult C. hongkongensis (2 year-old) were collected from Zhanjiang (Guangdong, China) in 2010. One hundred oysters of each species were acclimated in an aquarium tank supplied with sand-filtered seawater at ambient temperature (18 ± 1 °C) and salinity (25‰). After acclimation for a week, the two species of oysters were divided into two groups (control group and treatment group), respectively. Control groups were kept in filtered seawater with optimal salinity 25‰, while treatment groups were exposed to low osmotic filtered seawater with salinity of 8‰. In order to ensure the free exchange of seawater between the inside and outside of the shell, a part of the shell edge (about 10 mm long and 5 mm wide) of each specimen was chipped away. After 8 hours, six oysters from each group were randomly selected, and the gills were dissected for RNA extraction.

Small RNA library preparation and sequencing

Experimental protocols for the cDNA normalization sequences were performed according to the manufacturer’s technical instructions. Firstly, the total RNA was isolated from gill tissues with TRIzol reagent (Invitrogen), and the RNA of 15–40 nt was purified by the 15% Polyacrylamide gel. Equal amounts of the high-quality small RNA samples from six individuals of each group were then pooled for cDNA library preparation using the TruSeq Small RNA sample Preparation Kit (Illumina), respectively. Then these small RNAs were ligated sequentially to 5′ and 3′ adapters and used SuperScript II Reverse Transcriptase (Life Technologies) to synthesize the first-strand cDNA. The double-strand cDNA was performed followed by PCR amplification.

The purified PCR products were used for cluster generation by TruSeq PE Cluster Kit (Illumina) and then was sequenced on Illumina Hiseq 2000 following the manufacturer’s instructions. Raw reads obtained from Illumina sequencing have been deposited in the National Center for Biotechnology Information Short Reads Archive ( http://www.ncbi.nlm.nih.gov/srawebsite) under accession number SRP049540.

Data processing

The raw reads were subject to the program FastQC79 ( http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to assess the quality of sequencing data and trim low quality reads and adaptor sequence. The reads with high sequencing quality and ranging from 15–30nt in length were annotated by searching against the GenBank database80 and the Rfam database81. In addition, the clean reads of C. gigas were aligned to the Pacific oyster genome (GenBank accession: AFTI00000000.1, http://metazoa.ensembl.org/Crassostrea_gigas/Info/Index) while the reads of C. hongkongensis were aligned to the transcriptome of the Hongkong oyster48 using bowtie software82 with no more than 2 nt mismatches to filter the reads.

Analysis of conserved and novel miRNAs

Reads mapped to either non-miRNA in Rfam (such as rRNAs, tRNAs, snoRNAs etc.) or oyster mRNAs, were excluded for further analysis. The remaining reads were aligned against the miRBase 20.083 and the oyster genome for conserved and novel miRNA identification combined with stem-loop structure prediction. The clean reads were mapped to mature miRNA and hairpin sequences in miRBase 20.0 with complete matches to identify conserved miRNAs. Reads that did not match miRBase database were marked as unannotated. The unannotated data sets were aligned with the oyster genomic sequences to predict novel miRNAs, using miRDeep2 software84 with the prediction of the secondary structure.

Differential expression analysis of miRNAs

To compare miRNAs expression data between the control and treatment groups, read counts for each identified miRNA were normalized to the total number of reads in each given sample. Two methods, DESeq85 and edgeR86, were used to determine the differential expression among the experimental groups. These two methods are all have type-I error control87. EdgeR detected differentially expressed transcripts at lower sensitivity to DESeq, while DESeq returned more false-positive transcripts88. A common practice to boost the result accuracy is to use more than one method and then to combine their results. So only differentially expressed miRNAs identified by both methods were considered for further analysis. The fold change of miRNAs was calculated as the ratio of read counts in the treatment group to the read counts in the control group followed by transformation of log2. The miRNAs with the absolute value of log2 (fold change) ≥1.0, and q values <0.05 were considered as significantly differentially expressed.

Target gene prediction and functional annotation

The 3′UTR sequences of oyster protein-coding genes were retrieved based on the oyster genomic sequences and annotation information. The target genes of oyster miRNAs were predicted using the miRanda algorithm (score ≥160, free energy ≤−25 kcal/mol).

The target sequences were annotated using Blast2GO software89 for assigning GO terms to investigate their putative functions. The GO terms were plotted through Web Gene Ontology Annotation Plot (WEGO)90. Moreover, the target genes were also annotated using Kyoto Encyclopedia of Genes and Genomes91. To determine the possible overlapping of biological functions among these miRNAs significantly overrepresented GO terms and KEGG pathways were searched using the GOstat92 package and the GSEABase package93 with a P-value cutoff of 0.01.

Based on the results of miRanda algorithm, the differentially expressed mRNA from transcriptome data48 were compared with the target genes of differentially expressed miRNAs. The genes that were found to be differentially expressed in C. gigas adults in response to six different salinity treatments when compared to a control salinity of 30 Table S21 in Zhang et al.49 were considered as osmoregulatory candidates35. Additionally, we quantified the number of target genes of differentially expressed miRNAs that were included in osmoregulatory candidate genes.

qRT-PCR validation of miRNA expression

The miRNA were extracted and purified from gill tissues using High Pure miRNA Isolation Kit (Roche) according to the manufacturer’s instructions. The reverse transcription was carried out using on miScript II RT Kit (QIAGEN). The synthesis reaction was incubated for 60 min at 37 °C, and terminated by heating at 95 °C for 5 min to inactivate enzyme reaction.

The qRT-PCR was performed using the miScript SYBR Green PCR Kit (QIAGEN). The reactions were carried out in a total volume of 25 μl containing 2.5 μl of diluted cDNA, 2.5 μl of each primer, and 12.5 μl SYBR Green PCR Master Mix with the following cycling profile: 95 °C for 15 min for polymerase activation, followed by 45 cycles at 94 °C for 15s, 55 °C for 30 s, 70 °C for 30 s. Eight miRNA fragments were amplified using specific forward primers (Table 5) and universal reverse primers and 5S fragment were used as an internal control. Each sample was processed in triplicate and conducted with Roche Lightcycler 480 (Roche). All data was analyzed using 2 –ΔΔCt method.

Table 5

The primer sequences used for qRT-PCR in this study.

Primer	Sequence(5′-3′)
5S	TTGGATGGGTGACCGCCTG
miR-2353 (C.gigas)	TCTGTACTGCAGAATATCCAATATC
miR-1984	TGCCCTATCCGTCAGTCG
miR-92	GCCGGGACGAGTGCAATA
scaffold43364_10952	GTGCTCATTTGTCGAAACTGT
miR-183	TGGGAATGGCACTGGTAGAATTC
miR-184-3p	CCTTATCAGTTCTCCGTCCA
miR-3205	CACTCGGTGGACTGCTCG
miR-2353 (C.hongkongensis)	ACACTTCTCCAGTACAGAGGATCT

Conclusions

This is the first report to investigate the expression profiles of small RNAs in response to acute hypo-osmotic stress in oysters. A total of 202 and 87 miRNAs were identified from C. gigas and C. hongkongensis, respectively. Differential expression analysis of the miRNAs suggested that miRNAs played conspicuous roles in response to low salinity stress. Based on the GO annotation of target genes, miRNA mainly participate in the biological processes including microtubule-based process and cellular component movement, etc. Meanwhile the different expression patterns of miRNAs and the miRNA-mRNA interaction pairs indicated the differences in adaptation to hypo-osmotic stress between two oyster species. Further studies are required to understand the biological functions of miRNAs, and focus on the responsible miRNAs and Mollusca-specific miRNAs in response to osmotic stress. Additionally, increasing experimental sets would generate more conclusions.

Additional Information

How to cite this article: Zhao, X. et al. High throughput sequencing of small RNAs transcriptomes in two Crassostrea oysters identifies microRNAs involved in osmotic stress response. Sci. Rep. 6, 22687; doi: 10.1038/srep22687 (2016).