Characterization of microRNAs in Mud Crab Scylla paramamosain under Vibrio parahaemolyticus Infection

Two mixed RNA pools of 7 tissues (intestine, heart, liver, gill, brain, muscle and blood) from the V. parahaemolyticus infected mud crabs and the control groups, respectively, were subjected to high-throughput sequencing with the Illumina/Solexa platform. A total of 19,144,358 raw reads from the control group and 18,559,070 under the infected conditions, were obtained and uploaded to Gene Expression Omnibus (GEO: GSE39921↗). After removing some substandard reads, we analyzed the length distribution of the rest small RNAs (13–27 nt long) (control, 17,496,577, 91.90% of total, GEO: GSM981350↗; infected, 16,888,076, 91.51% of total, GEO: GSM981351↗) (Figure 1). The size distribution was similar between these two groups, and the most abundant small RNA was 22 nt long (Figure 1A). The small RNA sequences were compared with RNA databases (miRBase 19.0, http://www.mirbase.org↗, released at August 2012; Rfam, http://rfam.janelia.org↗; and Genbank, http://blast.ncbi.nlm.nih.gov↗). Reads from mRNA, rRNA, tRNA, snoRNA, and repeat sequences were removed. A quantity of 13,926,242 (^~80%) and 13,138,560 (^~78%) meaningful reads in the control and infected conditions, respectively, were remained for miRNA analysis (Figure 1 B, C).

To identify conserved miRNAs in mud crab, all small RNA sequences were Blastn searched against the currently known miRNAs and clustered on the basis of sequence similarity because some sequences probably came from the same precursor. To get a comprehensive knowledge of the conserved miRNAs in S. paramamosain, we compared mud crab miRNA transcriptome during V. parahaemolyticus infection and the control groups (Table S1) to known miRNAs in miRBase 19.0. Allowed one or two mismatches between sequences, after Blastn searches and further sequence analysis, total miRNAs in this two groups were belong to 421 miRNA families identified by miRBase (Table S2). Among the 421 families, 98 (23.3%) Spa-miRNA families were conserved in Arthropoda, 35 (8.3%) of which were conserved in Crustacea. Thus, 133 conserved miRNA families in mud crab S. paramamosain were either identical or very similar to existing miRNAs in miRBase, whereas the remaining 288 were known but non-conserved miRNAs in the crab miRNA transcriptome.

To validate the data acquired from sequencing, northern blotting analyses were conducted to determine the expression levels of 4 randomly-chosen significantly differentially expressed miRNAs (up-regulated: miR-266, miR-1778; down-regulated: miR-702-5p, miR-84a) with 5S rRNA gene of the mud crab as the internal reference. The results of the mixed RNA samples by northern blotting were well in accord with that of the sequencing by Solexa (Figure 2a), indicating the data presented here were in high quality. Furthermore, we measured the expression profiles of these four miRNAs by stem-loop quantitative real-time PCR assay in seven individual tissues. The data were subjected to analysis of one-way ANOVA followed by Student’s t-test (Figure 2c–f). miR-266 was expressed in all seven tissues tested in the treatment groups and the control groups with the highest expression level in the tissues of muscle, blood and liver. After V. parahaemolyticus infection, miR-266 increased its reads significantly in the muscle of the mud crab. Another up-regulated miR-1778, which had a very low reads in all tissues tested in the control groups, had a strong expression levels in blood and followed by that in muscle. As for two down-regulated miRNAs, miR-702-5p decreased its reads, most significantly, in muscle and liver; while miR-84a reduced its expression most in brain and gill after the V. parahaemolyticus infection in mud crabs. In summary, 4 differentially expressed miRNAs, all identified by Solexa sequencing, were validated by qRT-PCR assay, and the data showed a tissue specific characteristic only for miR-1778 (blood and muscle-specific).

To better understand the function of the above novel as well as conserved mud crab miRNAs, putative targets of these miRNAs were predicted using the described criteria and methods [24]. All miRNAs were used for target analysis by search against D. pulex EST sequences and the annotated cDNA sequences of mud crab S. paramamosain (unpublished data). All miRNAs have multiple targets, and some even have more than 2000. Similarly, some Unigenes have varied miRNA target sites. Top 3 targets by the order of the scores were obtained for each miRNAs, and subjected to GO analysis to investigate gene ontology. The results demonstrated that the targets of miRNAs in mud crab could be involved in biological processes, cellular components and molecular functions (Figure 3).

Moreover, 8small RNAs were identified in the crab miRNA expressions of the infection group and the control group from precursors by using Mireap, designed by Beijing Genomics Institute (BGI), met the criteria established by Allen et al. [25] and were considered as putative novel S. paramamosain miRNAs (Table 1). Since the whole genomic sequence of mud crab is unavailable, we rely on the genome sequences of crustacean Daphniapulex and our previous mud crab transcriptom database (unpublished data) as references to predict mud crab new miRNAs. The precursor sequences of the 8 novel potential miRNAs were verified by PCR amplification with S. paramamosain genomic DNA as templates (Figure S1), cloned and sequenced with the length varying from 76 to 100 nt (Table 1). The secondary structures were then predicted by Mfold to be a typical hairpin (Figure 4). Further, the minimum free energy (Mfe) value of the five miRNAs was estimated to range from -41.9 to -20.82 kcal/mol (Table 1), which were consistent with the results from previous researches [26,27].

The expression profiles of the 8 novel miRNAs during the pathogen infection groups and the control groups were assayed using qRT-PCR analysis and the signals were detected for all of them (Figure 2b). The stem-loop RT-PCR is a reliable method and a convenient way for detecting and identifying the novel miRNAs. In this study, we adopted this technique to validate and measure the expression of 8 novel potential miRNAs. All of these potential miRNAs were identified in S. paramamosain by Solexa sequencing. To aid in determination of S. paramamosain novel miRNA functions, we examined their expression in two libraries RNA pools. Of these 8 miRNAs, only Spa-novel-5 were significantly up-regulated, and the other 7 miRNAs were down-regulated, 6 of which were significantly reduced during the V. parahaemolyticus infection in S. paramamosain (Figure 2b). The data obtained here can form powerful evidence to support the existence of the novel miRNAs in S. paramamosain. The amplification plots and the melting plots of qRT-PCR analysis were provided as supporting information attached in Figure S2.

To identify the infection-responsive miRNAs from S. paramamosain, the number of normalized miRNA reads of the treatment groups and the control groups were compared. Based on the sequencing results, the differential expression of miRNAs greater than two-fold were chosen for calculation. Of the differentially expressed miRNAs, 288 (48.6%) were expressed in both the control and infected mud crabs, whereas 196 (24.4%) and 151 (27.1%) were specifically expressed in the control and infected groups, respectively. After normalization, the expression levels of 161 miRNAs (161/421, 38%) (106 up-regulated and 55 down-regulated) were found to be significantly different between these two samples. The scatter plot for differentially expressed miRNAs between the treatment groups and the control groups was showed in Figure 5.

Three distinct miRNA families were found to change dramatically in the infection groups: miR-1, miR-2 and let-7. The miR-1 family, including miR-1, miR-1a, miR-1b, miR-1c, and miR-1*, had the highest total copy number of 18,356,231. The miR-2 family was the second most highly represented with 8 members including miR-2a-f, miR-2, miR-2b*, and miR-2f*. The let-7 family had 13 members including let-7, let-7a–k, and let-7f* (For more information, please see). Table S1

The function of the differential expressed miRNA’s target is listed in Table S5. Actually, there are multiple targets for every miRNA. Meanwhile, for some genes, everyone has more than one miRNA binding sites. Part of some miRNAs targets would be immune related genes, such as miR-1345 targeting cryptocyanin and miR-135a* targeting kelch domain containing 2-like protein, Anoctamin-7 for miR-1224, calpain B for miR-1275, Immunoglubin I domain-containing protein for miR-412*, etc. Other predicted targets were annotated as transcriptor factor, such as transcriptase for miR-638, transcriptase-like protein for miR-1406 and transcription domain-associated protein for miR-198. Notably, several miRNAs targets would function in disease resistance, such as disease resistance protein (NBS-LRR) for two miRNAs miR328* and spa-novel-02. Due to the incomplete annotation of the transcriptomic data of S. paramamosain, the target of some miRNAs were function unkown, Moreover, The substantial part of the miRNAs identified during the infection in mud crab have no target genes at all. miRNA targets prediction and the function of the differential expressed miRNAs targets need further investigations.

Two conserved miRNAs (and respective homologoues), miR-965 (miR-965*, miR-965-3p) and miR-33 (miR-33a, miR-33b, miR-33-5p, miR-33a-5p, miR-33b-5p), were selected for further phylogenetic analysis (Figure 6). For miR-965 and its two homologoues, the seed regions, 1-9 and 11-19, were identical in mud crab and other crustaceans (D. pulex, shrimp Parhyalehawaiensis and mitten crab Spiroplasmaeriocheiris). However, at the 10^th position, a transition mutation of A to G was found in miR-965 of mud crab when compared to that of non-crustacean animals (Figure 6A). For miR-33 and its five homologoues, instead, a transition mutation of G to A at the 1^st position could be identified only in mud crab (Figure 6B). Therefore, nucleotide sequences of miRNAs of S. paramamosain not only show the similar evolutionary patterns within the crustaceans, but also possess mud crab-specific mutations.

The mud crabs used in this study were taken from a local crab farm (Niutianyang, Shantou, Guangdong, China). No specific permits were required for the described field studies, as the sampling locations were not privately owned or protected in any way. Furthermore these field studies did not involve endangered or protected species. The animals were processed according to "the Regulations for the Administration of Affairs Concerning Experimental Animals" established by the Guangdong Provincial Department of Science and Technology on the Use and Care of Animals.

Forty healthy mud crabs (100 ± 10 g in weight) were collected from Niutianyang area in Shantou, China. Crabs were acclimated for 3 days in 1 m³ tanks (ten crabs per tank) in conditions (salinity: 8‰; temperature: 28^oC) similar to those of culture ponds from which the crabs were obtained. Seawater was changed twice a day. Crabs were then fed with shellfishes for one week before sampling. V. parahaemolyticus isolated from diseased crabs (from Shantou crab farms) were cultured in 2216E medium at 28^oC. Twenty crabs were inoculated individually in the ventral hind legs with 100 µL of V. parahaemolyticus (2.4 × 10⁷ cells/mL), while the other twenty crabs were injected with 100 µL of saline (0.9%) as the control group. 4 hours after the injection of V. parahaemolyticus, crabs showed a frail situation of decreased appetite, reduced ingestion, paralyzed ambulatory legs, discolored carapace, and milky appendages. 20 crabs from each of the two groups were collected for RNA isolation 24 hours after the challenge. Haemocytes were collected from blood samples after mixing with anti-coagulant solution and centrifugation at 800 g at 4^oC for 20 min. The haemocytes and six other dissected tissues (intestine, heart, liver, gill, brain and muscle) were immediately frozen in liquid nitrogen and then stored at -80^oC prior to RNA extraction.

Total RNA was extracted using TRIzol® (Invitrogen, Carlsbad, CA) according to the manufacturer’s protocol. Genomic DNA was removed by RNase-free DNase I digestion to gain high quality RNA samples, which were confirmed through measuring the ratio of 260/280 nm and the RNA integrity number (RIN) value [42]. After quantification, equal amount of RNA samples (2 µg) from the 7 tissues were mixed before being used for small RNA library construction.

Small RNA fractions (18-30 nt) of the two mixed RNA samples under normal and infected conditions were enriched by polyacrylamide gel electrophoresis, and ligated with adapters to both 5'- and 3'-termini. Ligation products (62-75 nt) were reversely transcribed and those cDNAs were used for deep sequencing by an Illumina/Solexa G1 sequencer.

Sequences of 50 nt were obtained via the Illumina/Solexa sequencing platform of Beijing Genomics Institute (Shenzhen, China). After removing of adapters, low quality tags, and contaminants, clean sequences were gained for further study. By aligning these sequences with data in miRBase 18.0 website, Genbank website, and Rfam database, small RNAs were separated into six groups: miRNA, rRNA, tRNA, unknown sequences, snRNA/snoRNA and degradation fragments of mRNAs. A total of 19,144,358 raw reads from the control group and 18,559,070 under the infected conditions have been deposited in Gene Expression Omnibus (GEO; accession number GSE39921↗). After removing the substandard reads, 13,926,242 reads from the control group and 13,138,560 reads from the infected group have also been deposited with GEO accession number (control: GSM981350↗; infected: GSM981351↗).

Known miRNAs were elucidated by comparing genome sequences of Daphniapulex and our previous mud crab transcriptom database (unpublished data). To determine conserved miRNA, the filtered sequences were initially used to search the miRNA database, miRBase 19.0 with BLASTN. Only matching sequences (0-2 mismatches) were considered as conserved Spa-miRNA. Small RNA sequences were aligned to our previous mud crab transcriptom database to obtain precursor sequences in order to identify novel miRNAs. Contexts of perfectly matched hits were extracted as ±150 bp. The Mireap program developed by the Beijing Genome Institute (BGI) was used to analyze structural features of miRNA precursors to identify novel miRNA candidates. The resulting structures were retained as novel potential miRNAs only if they met the criteria described by Allen et al. (2005) [25]. The secondary structures of filtered pre-miRNA sequences were checked using Mfold [43]. In each case, only the structure with the lowest-energy was selected for manual inspection.

The sequence reads of the two libraries were normalized to 1 million by the total number of sRNA reads in each sample. The calculation of the p-value for comparison of the miRNA expression between the two libraries was based on previously established methods [44]. Specifically, the log₂ ratio formula was: log₂ ratio = log₂ (miRNA reads in infection treatment/miRNA reads in control).

The following p-value formulae were used:

Where N1 is the total number of reads in the sequencing library of the control, N2 is the total number of reads in the sequencing library of the infection treatment, x is the number of reads for an miRNA in the control library, and y is the number of reads for an miRNA in the infection treatment library. All calculations were performed on a BGI Bio-Cloud Computing platform (http://www.genomics.cn/en/navigation/show_navigation?nid=4143↗). Normalized miRNAs of <1 were filtered in both libraries.

The unknown sequences were utilized to forecast novel miRNA by Mireap software, which could be used to analyze miRNA precursor’s structural features and recognize novel miRNA candidates. The secondary structures of alternative pre-miRNA sequences were verified by Mfold software using a minimum free energy method.

The assay to amplify miRNA precursors was conducted as previously described [45]. Briefly, we extracted genomic DNA from the muscle tissue of the mud crabs according to the manufacture’s protocol and designed primers for 20 of 57 duplex-like miRNA: miRNA* pairs with Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA). PCR reactions were performed with the pfu enzyme. Corresponding fragments were amplified by PCR and the length of amplification products was detected on 3% agarose gels. Fragments between 60 and 100 nt in length were subcloned into pMD19-T vector (Takara Bio, Dalian, CA, China) for sequencing analysis.

Stem-loop RT-qPCR was used to confirm 8 novel potential miRNAs. The tissue distribution of 4 randomly-chosen differentially expressed miRNAs (miR-266, miR-1778, miR-702-5p and miR-84a) were also investigated by stem-loop RT-qPCR. All designed primers had been reported by Chen et al. [46] and listed in Table S6. cDNA was reverse-transcribed by PrimerScript ^®RT reagent Kit with gDNA Eraser (Perfect Real Time, Takara, Dalian, China) using a method previously described [47]. We performed quantitative real-time PCR using the SYBR^®premix Ex Taq^TM (Takara, Dalian, China) on an ABI Prism 7300 Sequence Detection System (Applied Biosystems, Foster City, CA). The reaction for each sample was performed in triplicate. Negative PCR products (no cDNA template) were prepared to detect any possible contamination. The specificity of primer amplicons was tested by analysis of a melting curve. The crab 5S rRNA was amplified as a reference gene for normalization using 2^-ΔΔT method [41,48,49]. All data were given in terms of relative miRNA expression as means±SE. The data were subjected to analysis of one-way ANOVA followed by Student’s t-test using SPSS 11.0 software, and the p-values smaller than 0.05 were considered statistically significant.

Northern blotting analyses were performed following the procedure as described in miRNA Northern blot user manual (Signosis, Inc., Sunnyvale, CA). Briefly, approximately 10µg of total RNA was loaded onto a 15% polyacrylamide urea denature gel in 0.5% TBE and run at 70V for about 45 min. After transferred in 0.5% TBE onto a nylon membrane at 70V for 60 min, RNA samples were fixed on the membrane by UV cross-linking. The membranes were hybridized with 10µl of biotin-labeled miRNA-1778, miRNA-84a, miRNA-702-5p and miRNA-266 probe overnight respectively. The membranes were then rinsed with 1 X hybridization washing buffer, and subsequently hybridized with 8 µl of the amplifier for two hours at 42 ^oC. After washed with 1x hybridization washing buffer for 30 min at 42 ^oC, the membranes were blocked with blocking buffer at room temperature for 30 min, then detected with Streptavidin-HRP. After washing three times with 1x Detection washing buffer, the membranes were overlaid with substrate A and substrate B for 5 min. The image was acquired using a FluorChem FC2 imager (Alpha Innotech Corp, USA). The RNA load was normalized by mud crab 5S rRNA.

Previous studies have proved that all miRNAs fulfilling the function of post-transcriptional gene regulation by binding to the target mRNA sequences in one or more perfect or near-perfect complementary site (s), bring convenience to predict animal miRNA targets simply using gene-homology search. To explore the potential function of miRNAs with significantly differential expression during the treatment with V. parahaemolyticus challenge in mud crab, we used the newly identified 161 differentially expressed miRNAs (known and novel miRNAs) to search their putative targets against annotated mud crab transcriptome data and published sequences from other invertebrates with RNAhybrid [24].