Differential expression of microRNAs in hemocytes from white shrimp Litopenaeus vannamei under copper stress

Fish and Shellfish Immunology 74 (2018) 152–161

Contents lists available at ScienceDirect

Fish and Shellfish Immunology journal homepage: www.elsevier.com/locate/fsi

Full length article

Differential expression of microRNAs in hemocytes from white shrimp Litopenaeus vannamei under copper stress

T

Hui Guoa, Zhi-cheng Lua, Xiao-wen Zhua, Chun-hua Zhua, Cheng-gui Wanga, Yu-chun Shena,∗, Wei Wangb,∗∗ a

Key Laboratory of Marine Ecology and Aquaculture Environment of Zhanjiang, College of Fisheries, Guangdong Ocean University, Zhanjiang 524025, People's Republic of China b College of Fisheries, Guangdong Ocean University, Zhanjiang 524025, People's Republic of China

A R T I C L E I N F O

A B S T R A C T

Keywords: MicroRNA Litopenaeus vannamei Cu Hemocytes

MicroRNAs (miRNAs) are small noncoding RNAs that regulate diverse cellular processes, including organismal stress response, through posttranscriptional repression of gene transcripts. They are known to have antiviral functions in aquatic crustacean species, but little is known about the role of miRNAs against environmental stress caused by Cu, a common chemical contaminant in aquatic environment. We performed small RNA sequencing to characterize the differentially expressed microRNAs in Cu exposed shrimp. A total of 4524 known miRNAs and 73 novel miRNAs were significantly (P < .05) differentially expressed after Cu exposure. The peak size of miRNAs was 22 nt. Among them, 218 miRNAs were conserved across 115 species. The validation of 12 miRNAs by stem-loop quantitative RT-PCR were found to be coherent with the expression profile of deep sequencing data as evaluated with Pearson's correlation coefficient (r = 0.707). Target genes of these differentially expressed miRNAs related to immune defense, apoptosis, and xenobiotics metabolism also showed significant changes in expression under Cu stress. The present study provides the first characterization of L. vannamei miRNAs and some target genes expression in response to Cu stress, and the findings support the hypothesis that certain miRNAs along with their target genes might be essential in the intricate adaptive response regulation networks. Our current study will provide valuable information to take an insight into molecular mechanism of L. vannamei against environmental stress.

1. Introduction

been reported to decrease the total hemocyte count (THC), respiratory burst activity and phenoloxidase activity of hemocytes in L. vannamei [7]. In our previous study, exposure of L. vannamei to Cu induced hemocyte apoptosis and the expression of antioxidant biomarker genes, apoptosis-related genes and a specific biomarker gene of heavy metal pollution [8]. However, the molecular mechanism of Cu stress in shrimp remains largely unclear. Biological systems use a variety of mechanisms to maintain their functions in the face of environmental and genetic perturbations [9]. MicroRNAs (miRNAs), a class of small noncoding RNAs with 20–24 nt in length, can regulate gene expression at the post-transcriptional level [10] and play crucial regulatory roles in a large variety of biological processes, such as development, cell differentiation, apoptosis, oncogenesis, immune and stress response in various organisms [11–14]. In animals, miRNAs regulate gene expression through imperfect sequencespecific binding to the 3′-untranslated regions (3′UTR) of the target mRNAs and usually cause translational repression [15]. A growing

Litopenaeus vannamei is the most important farmed shrimp species and is economically valuable. With the development of industrialization and intensive aquaculture, aquatic pollution has become a severe and growing problem. Among various environmental pollutants, heavy metals have been of great concern because of their inherent toxicity, non-degradability and persistence [1]. In aquaculture practices, copper sulfate is commonly used to eradicate filamentous algae and phytoplankton in shrimp farms. The excess application of copper sulfate in pond management may deteriorate the aquatic environment due to Cu accumulation in the pond sediments [2,3]. In crustaceans, hemocytes in the circulating hemolymph have vital roles in immune responses against pathogens, such as phagocytosis, coagulation, encapsulation, nodulation, production of antimicrobial peptides (AMPs) and proPO activity [4,5]. Consequently, toxic effects on hemocytes potentially affect the survival of these animals [6]. Cu has ∗

Corresponding author. Corresponding author. E-mail addresses: [email protected] (Y.-c. Shen), [email protected] (W. Wang).

∗∗

https://doi.org/10.1016/j.fsi.2017.12.053 Received 29 September 2017; Received in revised form 19 December 2017; Accepted 28 December 2017 Available online 02 January 2018 1050-4648/ © 2017 Published by Elsevier Ltd.

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

into an individual eppendorf tube held on ice, and then samples of hemolymph from three replicates were pooled for RNA extraction. The total RNAs of hemocytes were extracted using miRNA isolation kit (Ambion, USA) according to the manufacturer's protocol. The integrity and concentration were assessed using the Agilent Bioanalyzer 2100 System (Agilent Technologies, USA). The average RIN (RNA integrity number) value of samples was 8.9. Subsequently, the small RNAs ranging from 18 to 30 nt were separated by denaturing 15% polyacrylamide gel. After recovery by ethanol precipitation, the small RNAs were ligated sequentially to a pair of RNA adapters both in 5′ and 3′ ends by combing the reagents and incubating on the preheated thermal cycler according to the manufacturer's instructions. The reverse transcription and polymerase chain reaction (PCR) amplification were preformed after the ligation. The last products were sequenced using the Illumina Genome Analyzer (Illumina, San Diego, CA, USA) following the manufacturer's protocol and instructions.

number of miRNAs were discovered and studied in variety of organisms since the first miRNA lin-4 has been discovered in Caenorhabditis elegans. To date, over 28 thousand miRNAs have been discovered across 223 species (miRBase, release 21.0; June 2014). Initially, miRNAs were cloned and identified by a traditional method in M. japonicus [16]. However, PCR cloning followed by traditional sequencing remains very labor- and cost-intensive with a limited dynamic range to detect and define relative miRNA expression [17]. Recently, next-generation sequencing (NGS) technologies, including the Illumina Genome Analyzer (GA), Applied Biosystems SOLiD System, and 454 Life Sciences (Roche) FLX instruments, have emerged as well-established approaches for miRNA profiling, which can rapidly produce millions of sequence reads simultaneously with lower cost than Sanger sequencing and allows for the identification of low abundance miRNAs and genome-wide discovery of novel tissue-specific miRNAs with high specificity and sensitivity [18–20]. Deep-sequencing technologies have facilitated a sharp rise in the rate of novel microRNA discovery [21]. In recent years, miRNA studies have been performed using NGS technologies in crustaceans with limited genomic information [22–24]. Although hundreds of miRNAs have been identified, only a small number of shrimp miRNAs have been discovered and functionally identified in Marsupenaeus japonicus [22], Penaeus monodon [25] and L. vannamei [23] that responded to Vibrio alginolyticus or white spot syndrome virus infection. The investigation of miRNAs that responds to copper stress has not yet been performed. The objective of this study was to identify and characterize the differentially expressed miRNAs in shrimp L. vannamei in response to copper stress. The results will extend the knowledge of crustacean miRNA regulation and provide information for further research on shrimp response against environmental stress.

2.4. MiRNA identification and expression analysis The raw reads obtained from HiSeq were processed by removing low quality reads. After removal of the reads with proportion of N (N means unable to determine the base information) is greater than 10, 5′ primer contaminants and without the insert fragments reads, without 3′ primer reads, poly A reads and reads shorter than 18 nt, the remaining reads were regarded as clean reads and were blasted against the Rfam database (http://www.sanger. ac.uk/software/Rfam) and the GenBank noncoding RNA database (http://blast.ncbi.nlm.nih.gov/) to identify rRNA, tRNA, snRNA, snoRNA, mRNA, and repeat sequences. Those sequences were discarded assuming that they represent degradation and other undesired products. As there was no genome data available for L.vannamei, the remaining sequences were searched against the miRNAs from all the animals in miRBase 21.0 (http://www.mirbase. org/) using BLAST to identify known miRNAs. The steps are as following: (1) considering the difference among species, align clean data to the miRNA precursor/mature miRNA of all animals in miRBase allowing two mismatches and free gaps; (2) choose the highest expression miRNA for each mature miRNA family which is regarded as a temporary miRNA database; (3) align clean data to the above temporary miRNA database and the expression of miRNA is generated by summing the count of reads which can align to the temporary miRNA database within two mismatches; (4) predict the precursor of the identified miRNA, and those that cannot fold stem-loop hairpin structure will be regarded as pseudo-miRNA. The feasibility of the result can be greatly improved by this verification. The novel miRNA candidates were predicted using miRdeep software [27] with L. vannamei transcriptome deep sequencing data [28] based on hairpin-like secondary structure pattern.

2. Materials and methods 2.1. Animals The experimental shrimp L. vannamei (4.94 ± 0.50 g) were obtained from a commercial farm in Zhanjiang (Guangdong, China), and acclimated for 2 weeks prior to experiment in cycling-filtered plastic tanks with aerated seawater at 26 ± 2 °C and a salinity of 5‰. During the acclimation period, shrimp were fed twice daily with shrimp diet (40% protein, 5.0% fat, 5.0% fiber and 16% ash, supplied by a commercial diet, China) until 24 h before the experimental treatments began. Only shrimp apparently healthy and in the intermolt stage were used for the study. The molt stage was determined by examining the uropoda in which partial retraction of the epidermis could be distinguished [26]. 2.2. Cu exposure

2.5. Differential expression analysis

Exposure experiment was conducted in triplicate with thirty shrimps in plastic tank with 180 L water (26 ± 2 °C, pH 7.9–8.0, and salinity 5‰). Each tank was aerated continuously using an aeration stone. Normal and experimental concentrations were 0 and 5 mg L−1 according to previous study [7]. Experimental Cu concentrations were prepared by adding copper sulfate to 5‰ seawater until the desired concentration was attained.

To find out the differentially expressed miRNAs between two libraries, the expression of miRNA in 0 h and 3 h libraries were normalized to obtain the expression of transcripts per million (TPM) [29] using the following formula: Normalized expression = (Actual miRNA count/Total count of clean reads)*1,000,000. Then, the fold-change (Fold change = log2(3 h/0 h) and P-value [30] were calculated from the normalized expression. Identification of differentially expressed miRNAs between control and exposed samples was performed by pairwise comparisons across 0 h and 3 h samples using the DESeq R package (http://bioinfo.au.tsinghua.edu.cn/software/degseq/), and a corrected P-value < .05 and a threshold absolute |log2(3 h/0 h)|≥1 were set as the threshold for significantly differential expression [31,32].

2.3. RNA extraction, library construction and sequencing At the beginning (0 h) and after 3 h of Cu exposure, three shrimp were randomly sampled from each tank (n = 9). Hemolymph (200 μL) was withdrawn from the ventral sinus of each shrimp by a 1 mL sterile syringe (25 gauge) containing an equal volume of ice-cold anticoagulant (glucose 20.5 g L−1, sodium citrate 8 g L−1, sodium chloride 4.2 g L−1, pH 7.5). The hemolymph from each shrimp was transferred 153

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

18 nt, a total of 25,686,384 (89.76% of the raw reads) and 21,304,558 (87.77% of the raw reads) cleans reads, containing 845,582 (0 h) and 1,467,075 (3 h) unique sequences, respectively, were generated from the two samples. All the clean reads were annotated using GenBank, Rfam, miRbase 21.0 and our L. vannamei transcriptome databases, and classified into six categories including, rRNA, tRNA, snRNA, snoRNA, miRNA and those without any annotation. The count and proportion of the different categories of small RNAs were showed in Table 2. Among the different types of small RNA, only a small proportion of small RNAs were annotated (11.54% and 9.81% for 0 h and 3 h, respectively). The length analysis of the small RNAs showed that the two libraries shared a similar size distribution, which both varied from 20 to 23 nt and exhibited a unimodal distribution pattern, with 22 nt being the most abundant class, followed by the 21-, 20-, and 23-nt classes (Fig. 1).

2.6. Target gene prediction and enrichment analysis The target genes of differentially expressed miRNAs were predicted with miRanda (http://www.microrna.org) [31] and PITA (https:// genie.weizmann.ac.il/pubs/mir07/mir07_dyn_data.html) [36] computational target prediction algorithms. miRanda was used to match the entire miRNA sequences, and the parameters were set as free energy < -20 kcal/mol and score > 50 [22]. PITA was used to search for miRNA seed matches based on target-site accessibility and free energy. L. vannamei transcriptome deep sequencing data were used as the database for the target gene prediction [28]. The final results predicted by the two algorithms were combined and the overlaps were calculated. Enrichment analysis of these predicted target genes were analyzed by Gene Ontology (GO) (http://www.geneontology.org/) and KEGG pathway (http://www.genome.jp/kegg/). The Bonferroni Correction was carried out to correct P-values, only GO terms with P ≤ .05 and pathways with FDR≤0.05 were identified as highly represented.

3.2. Identification of known and novel miRNAs in L. vannamei To characterize the shrimp miRNA homologs, the small sequences were blasted against metazoan miRNA precursors and mature miRNAs in the latest miRBase v21.0 (http://www.mirbase.org/). A total of 6014 and 5823 sequences were identified to be similar to known miRNAs from other species that had previously been deposited in miRbase, in 0 h and 3 h libraries, respectively. To date, although the miRBase v21.0 has included miRNAs from 223 species, miRNA members from nonmodel animals have not been fully identified. To analyze the conservation of shrimp miRNAs, we compared them to those identified in animals, including worm (Caenorhabditis_elegans), mouse (Mus musculus), chicken (Gallus gallus), fruit fly (Drosophila_grimshawi), owl limpet (Lottia gigantea), water flea (Daphnia pulex), human (Homo sapiens) and etc. Among these known miRNAs, 218 conserved unique mature miRNAs were identified, which have high confidence miRNA orthologs in 115 species, with miR-9, miR-1 and let-7 are highly conserved from worm to human (Supplementary File Table S1). The unannotated small RNAs were matched with L. vannamei transcriptome database to predict potentially novel miRNA using mirdeep software. A total of 62 and 121 putative novel miRNAs were obtained from 0 h to 3 h libraries, respectively.

2.7. Validation of miRNAs and target genes To evaluate the accuracy of expression profiling for the miRNAs and their target genes, the obtained results were validated by quantitative RT-PCR (qRT-PCR) analysis. For miRNA validation, stem-loop qRT-PCR was conducted. Total RNA was extracted from the hemocytes of shrimp at 0 h or 3 h using miRNA isolation kit as mentioned above. Experiments were performed in triplicates, and at least nine shrimp were analyzed after quantified using NanoDrop 2000 spectrophotometer, the miRNAs were subjected to reverse transcription by TaqMan microRNA Assays kit (Applied Biosystems, USA) as mentioned in the previous report [33]. The reverse transcription primers and realtime PCR primers of selected miRNAs and U6 were designed based on stem-loop miRNA primer design principle, and the sequences were showed in Table 1. A total of 8 selected target genes, targeted by miR107, miR-7, miR-92b, miR-750, miR-184 and miR-128 were used for validation by qRT-PCR according to our previous report [34]. Briefly, the qRT-PCR of miRNAs and mRNAs were carried out in an Applied Biosystems 7500 (Applied Biosystems, Life Technologies, USA) using SYBR Green Real-time PCR Master Mix (TOYOBO, Japan) following the manufacturer's recommendations, and the data were calculated according to the 2−ΔΔCt method. All reactions were performed in triplicate. U6 and β-actin were used as internal controls to normalize the experimental data of miRNAs and mRNAs qRT-PCR, respectively.

3.3. Differentially expressed miRNAs induced by Cu exposure To identify miRNAs involved in Cu exposure, the expression profiles of the known and novel miRNAs were examined at 0 h and 3 h. The results revealed that the expression levels of 4524 known miRNAs and 73 novel miRNAs were significantly (p < .05) changed by Cu exposure. A total of 4293 known miRNAs were up-regulated and 231 known miRNAs were down-regulated, while 44 novel miRNAs were upregulated and 29 novel miRNAs were down-regulated. The miR-1175a3p and novel-miR-46 were the most significantly (P < .05) up-regulated, while miR-228 and novel-miR-8 were the most significantly (P < .05) down-regulated in known and novel miRNAs, respectively (Supplementary File Table S2). Overall, the most differentially expressed miRNAs were up-regulated after the Cu exposure.

2.8. Pearson correlation coefficient Correlation is a technique for investigating the relationship between two quantitative, continuous variables. A correlation relationship of miRNA deep sequencing counts and stem-loop qRT-PCR was assessed with Pearson correlation coefficient (r), a measure of the strength of the association between the two variables. A scatter plot of the variables was drawn to check for linearity between two variables. The nearer the scatter of points is to a straight line, the higher the strength of association between the variables. The Pearson's correlation coefficient r may take any value from -1 to +1 [35].

3.4. Validation of miRNAs and target genes

3. Results

To validate the credibility of HiSeq sequencing and bioinformatics analysis results, 12 differentially expressed miRNAs including 3 randomly selected novel miRNAs and 9 known miRNAs were subjected to stem-loop qRT-PCR validation (Fig. 2A). The results suggested that most miRNAs identified in this study were credible. Predicted target genes of some miRNAs were selected to validate the expression pattern using qRT-PCR (Fig. 2B). The results showed that the expression level of immune related genes, such as α-2M, C type lectin (CTL) and toll protein, detoxification related genes such as glutathione S-transferase (GST) and cytochrome p450 (CYP450), heat shock proteins (Hsp60), apoptotic related proteins such as inhibitor of apoptosis protein (IAP)

3.1. Illumina sequencing of miRNAs To identify miRNAs involved in the response of L. vannamei to Cu exposure, two separate libraries of small RNA were constructed, with the mixed pools of hemocytes from shrimp after exposure for 0 h and 3 h, respectively. The two libraries were subjected to Illumina sequencing, and a total of 28,879,125 (0 h) and 24,495,187 (3 h) raw sequences were obtained, respectively. After removal of low quality reads, adapter sequences, poly A sequences, and sequences smaller than 154

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

Table 1 Primers used for qRT-PCR validation. ID

Primer

novel-miR-82

F: GCGCGCGTGAGTGAGTGAG R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCACTCA F: CGCGCGCGGCCGACGA R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACATTTCA F: CTCACGCGGCCGACGATGA R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGGTATT F: CGCGTTGTTCGTTCGGCTCG R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTTGACG F: GCGCGAGCAGCATTGTACAGG R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGATAGC F: GCGCGTGAAAGACATGGGTAGT R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACAATCTC F: GCGCGAATTGCACTAGTCCCG R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGCAGGC F: CGCGCCAGATCTAACTCTTCCAT R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCGTCAT F: CAGTCGGAGAAAGACATGGGTAG R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACATCTCA F: CGACGTGTGAGGTAGTAGGTTGT R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACAACTAT F: ATTCGTGTGGCAGTGTGGTTAGC R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCAACCA F: GCTTGAGTGGACGGAGAACTGAT R: ATCCAGTGCAGGGTCCGAGG RT:GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGCCCTT F: CTCGCTTCGGCAGCACA R: CTCGCTTCGGCAGCACA F: GCCCATCTACGAGGGATA R: GGTGGTCGTGAAGGTGTAA F: GCACGTAATCAAGATCCG R: CCCATCTCATTAGCACAAAC F: ATGCCAACGGTTTATTCGGAC R: GAGGTACAAAGCCCGCAACAC F: GCAGTTTGCACAGACGTGTT R: CCATGCTGGGAAAACGATTA F: AAGATAACGCAGAGCAAGG R: TCGTAGGTGACGGTAAAGA F: ATTGTCCGCAAGGCTATC R: ATCTCCAGACGCTTCCAT F: GAGATGACGTGAAATGCTTTGACT R: TGCCGGACTGATGGTCTTG F: CGAATCCCCACATCCACG R: GGCGGCTGATACACCACC F: TCGACCATCCCTTTTACACC R: TTGCCTGGAAGGTCTGATTC

novel-miR-48

novel-miR-47

miR-375

miR-107-3p

miR-7

miR-92b

miR-750-3p

miR-71-5p

let-7a-5p

miR-34-5p

miR-184

U6 β-actin α-2M

CTL CYP450 GST Hsp60 IAP p53 toll

and p53 were significantly changed by Cu exposure. Examination of the relationship between HiSeq deep sequencing and qPCR expression levels revealed a close correlation between the expression changes (fold difference) measured by each method (regression P < .05; Pearson's correlation coefficient r = 0.707) (Fig. 3).

Table 2 Summary of all small RNAs from two libraries of L. vannamei. Types

rRNA snoRNA tRNA miRNA snRNA unann Total

0h

3h

Number

Percentage (%)

Number

Percentage (%)

62485 797 28197 4208 1925 747970 845582

7.39 0.09 3.33 0.50 0.23 88.46 100

103282 941 28374 4847 6583 1323084 1467075

7.04 0.06 1.93 0.33 0.45 90.19 100

3.5. Target gene prediction and functional analysis of differentially expressed miRNAs The putative target genes of differentially expressed miRNAs were predicted using miRanda and PITA. A total of 21348 and 12579 target genes were predicted respectively for the 4524 differentially expressed known miRNAs and 73 differentially expressed novel miRNAs by the two prediction programs. In this study, many miRNAs have more than 155

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

Fig. 3. Comparison of relative fold change obtained from deep sequencing and qPCR analysis of 12 selected miRNAs. Each point represents the expression 3 h relative to 0 h, as measured by deep sequencing and qPCR. The regression indicates a close correlation between methods (linear regression of log2 fold differences: P < .05; r = 0.707).

Fig. 1. Length distribution of miRNAs in 0 h and 3 h libraries from L. vannamei.

one predicted target gene. Apoptosis related genes such as p53 protein, caspase, immune related genes such as CTLs, Cactus, Toll-protein, Spaetzle, stress response related genes such as Hsp60, Hsp70, Hsp21, detoxification related genes CYP450 and GST were targeted by some miRNAs (Supplementary File Table S3). The predicted target genes, which were stress response-related, of the identified L. vannamei miRNAs were categorized into 9 groups, such as ProPO system, oxidative stress, detoxification, proteinase/proteinase inhibitors, pattern recognition proteins, heat shock proteins, apoptotic related proteins, antimicrobial peptides and other immune molecules (Supplementary File Table S4). The predicted interactions of some miRNAs and target immune genes in Cu-exposed shimp are shown in Fig. 4. To explore the potential pathways regulated by these miRNAs, the target genes of the differentially expressed miRNAs were subjected to GO analysis and KEGG pathway analysis. According to GO classifications, 21348 and 12579 target genes, which were predicted for differentially expressed known and novel miRNAs, were clustered into 56 and 52 GO terms, respectively. The top ten terms according to Biological Process (BP), Cellular Component (CC) and Molecular Function (MF) classifications are illustrated in Fig. 5. Similar to GO analysis, KEGG pathway analysis was also based on target genes of known and novel miRNAs with differential expression. The results showed that 21348 target genes of known miRNAs were grouped into 307 pathways while 12579 target genes of novel miRNAs were grouped into 308 pathways. As we mapped the top 20 of enriched pathways in descending order of the rich-factor, we found the amoebiasis pathway (ko05146) and vibrio cholerae infection pathway (ko05110) seemed to be some important pathways under Cu stress (Fig. 6 and Fig. 7).

4. Discussion Living animals are constantly faced with various environmental stress that challenge normal life. In recent years, many studies in comparative biochemistry field are now on the verge of dissecting the role of miRNAs in adaptation to environmental stress [37]. It has been suggested that miRNAs serve as key players in a robust adaptive response against stress in animals through their capacity to fine-tune gene expression [38]. The role of miRNAs in environmental stress response has been studied in corals [39], Caenorhabditis elegans [38], Littorina littorea [40], Apostichopus japonicus [41] and other mammals. However, no comprehensive miRNA study has been performed on white shrimp under environmental stress. In this investigation, we revealed the

Fig. 2. Expression levels of selected miRNAs (A) and target genes (B) in the hemocytes of L.vannmei under Cu stress. These data are expressed as the mean ± SD relative to the control. *p < .05 vs. control; **p < .01 vs. control.

156

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

Fig. 4. The interactions between miRNAs and target immune genes of L.vannamei. Note: Triangles indicate miRNAs, blue circles indicate target genes. Triangle color represents different fold-change, ranging from pink (-0.33) to dark red (3.84). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

exposure for 0 h and 3 h to Cu stress. The results further confirmed that the de novo transcriptome has aided in the prediction of novel miRNA in crustacean [45] and reducing the possibility of over prediction [51]. To identify the miRNAs involved in the host response to the Cu exposure, the miRNA expression levels at 3 h were compared with the expression levels in controls. The results revealed that a total of 4524 known miRNAs and 73 novel miRNAs were significantly differentially expressed implying they are involved in the stress response. Some miRNAs belong to the family of miR-1175 and miR-750 were the most significantly up-regulated in the present study. This result was consistent with the previous study in L. vannamei after WSSV infection, suggesting these miRNAs were mainly related to cell function and signal transduction [52]. 3 novel miRNAs and 9 known miRNAs identified in our study were validated by stem-loop qRT-PCR, demonstrating a high level of confidence for the obtained data. The relationship analysis between deep sequencing and qRT-PCR revealed the qPCR results of 12 miRNAs agreed with the deep sequencing results (Pearson's correlation coefficient r = 0.707). The Pearson's correlation coefficient (r) was close to that of the previous studies [45]. The consistence further confirmed that the miRNAs identified by HiSeq and bioinformatics should be credible. The target genes of differentially expressed miRNAs were predicted using bioinformatic methods to characterize miRNAs function in shrimp stress response against Cu exposure. In the present study, the highestexpressed miRNAs were miR-1175-3p and novel-miR-46; and the lowest-expressed miRNAs were miR-228 and novel-miR-8. The prophenoloxidase (proPO) was a predicted target gene of miR-1175-3p. proPO has been extensively investigated in crustaceans due to its crucial role in non-self recognition as well as several subsequent immunodefense processes that are particularly essential to resist microbial infections [53,54]. Spz3, a predicted target gene of novel-miR-46, is a signaling ligand in innate immune response [55]. Relish is a predicted target gene of miR-228 and has been reported to be involved in immune-deficiency (IMD) signal pathway [56]. Kazal-type serine proteinase inhibitor, a predicted target gene of novel-miR-8, is an important group of serine protease inhibitors (SPIs) that can modulate several immune processes [57]. In addition to the highest- and lowest-expressed miRNAs, predicted target genes of some other significantly (P < .05) differentially expressed miRNAs were reported to participate

repertoire of miRNAs in L. vannamei under Cu stress through RNA sequencing and bioinformatics analysis. A total of 4524 known miRNAs, 73 novel miRNAs and 218 conserved miRNAs were identified after Cu exposure. In the present study, the peak size of miRNAs in L. vannmei was 22 nt. This result is consistent with previous studies in L. vannamei [23], F. chinensis [42] and M. japonicus [43] indicating that 22 nt might be a typical size of miRNAs in these crustacean species [42]. The miRNAs identified in this study can serve as the foundation for further research on shrimp miRNA, and will help to fill the knowledge gap in miRNA function in stress response of shrimp. The number of known miRNAs, conserved miRNAs and novel miRNAs identified in this study was much higher than that of WSSV [43] and Vibrio alginolyticusinfected M. japonicus [22], WSSV-infected P. monodon [25] and L. vannamei [23]. This might be due to the different miRBase version used in the studies. The miRBase database from version 15.0 to 19.0 was used in the previous studies mentioned above, whereas an update of the miRBase database (miRBase 21.0) was applied in the present study. The increasing number of miRNA sequences in the latest miRBase provide more comprehensive references for miRNA annotation, discovery and prediction. Many reported miRNAs are phylogenetically conserved [44]. All of the 5 miRNAs of M. japonicus in the miRBase database (http://www. mirbase.org/cgi-bin/mirna_summary.pl?org=mja) were matched with the miRNAs of L. vannamei identified in the present study. The result is in good agreement with the previous study that the 5 miRNAs were all identified in Fenneropenaeus chinensis [42]. Homologs of known miRNAs reported in the shrimp M. japonicus [43], P. monodon [25], M. rosenbergii [45] and L. vannamei [23], such as miR-750, miR-2001, miR2765, miR-2001, miR-1, miR-7, miR-34, miR-133, miR-317 and miR965, have also been identified in the present study, implying that these miRNAs might be highly conserved in shrimp. Some miRNAs such as miR-9 [25], miR-21 [46], miR-124 [47], miR-142 [48], miR-138 [49] which have already been described in vertebrates were found in the present study, suggesting these miRNAs may be conserved from vertebrates to invertebrates. The emergence of novel miRNA sequences may provide opportunities for combinatorial miRNA regulation that contribute to a greater potential for genetic robustness and morphological or functional complexity [39,50]. In the current study, a total of 62 and 121 novel miRNA candidates were predicted, respectively after 157

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

Fig. 5. Top ten enriched GO terms for the target genes of differentially expressed miRNAs in L. vannamei under Cu exposure.

The predicted target genes and the functions of their products provide useful clues for research on some essential biological processes in shrimp. According to GO analysis, the three most highly enriched GO terms in the BP, CC and MF categories were cellular process, cell and binding, respectively, which is consistent with the previous studies in liver and gill of Larimichthys crocea [58] and in hemocytes of Crassostrea gigas [59]. Some GO terms associated with immune and stress response, (such as MyD88-dependent toll-like receptor, regulation of adaptive immune response and response to stimulus), and repair (such as DNA repair, Oxidative demethylation and negative regulation of DNA repair) are also enriched in the target genes of differentially expressed miRNAs. These results were similar with the previous research in in Acropora digitifera under thermal stress [39]. Intriguingly, amoebiasis (ko05146) and vibrio cholerae infection (ko05110) were the most enriched

in shrimp immunity, oxidative stress, detoxification and apoptosis (Fig. 4). Some genes were predicted to be regulated by multiple miRNAs, and likewise, a single miRNA might target more than one genes. The results were consistent with the research of P. monodon challenged with WSSV infection, in which the immune-related target genes, apoptosis-related proteins, oxidative stress proteins and other immune molecules were also predicted [25]. Validation of several immune-related, apoptosis-related, or detoxification related target genes of some miRNAs showed that the expression levels of these genes were significantly changed by Cu exposure implying that these target genes might be involved in the shrimp defense mechanisms against environmental stress. However, further investigation is required to explore the potential function and molecular mechanism for these miRNAs and their target genes. 158

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

Fig. 6. The top 20 KEGG pathway analysis of known miRNAs with differential expression. Note: The horizontal coordinates is rich factor. The bigger rich factor the higher enrichment.

against Cu exposure. The predicted targets and their functional roles in stress response reported here indicate areas of interest that should be further investigated using independent methods. In conclusion, NGS technology and bioinformatics analysis was applied to identify miRNAs from hemocytes of L.vannamei under Cu stress. The results showed that 218 conserved miRNAs were identified, and the expression levels of 4524 known miRNAs and 73 novel miRNAs were significantly (p < .05) changed by Cu exposure. The deep sequencing results were further verified by a real-time RT-PCR assay. There was a consistency between the deep sequencing and real-time RTPCR method. Function annotation of the predicted target genes of the differentially expressed miRNAs revealed a broad range of pathways and processes related to immune defense, apoptosis, or xenobiotics metabolism were significantly enriched after Cu exposure. These findings will provide valuable information for further investigation on the mechanism of miRNA-mediated post-transcriptional regulation in the adaptive response of L. vannamei against environmental stress.

pathways in target genes of both known and novel miRNAs according to the KEGG analysis. In these two enriched pathways, the most enriched target gene was the MUC2/Mucin, an important component of innate defense in a number of enteric infections, including many parasitic infections [60]. It has been reported that during pathological conditions such as inflammatory bowel disease (IBD), MUC2 biosynthesis and secretion were accelerated [61], and was considered as a first step towards characterizing the innate immune response in different species [62]. The two significantly enriched pathways might indicate a potential defensive strategy of shrimp against Cu stress, since as a species living in water, shrimp may always face challenges of bacterial infection, especially under Cu stress, which could depress the immune ability and increase the susceptibility to bacteria infection [7]. Some immune related pathways such as JAK-STAT signaling pathway (ko04630), PPAR signaling pathway (ko03320), NOD-like receptor signaling pathway (ko04621), NF-κB signaling pathway (ko04064) and etc, drug and xenobiotics metabolism related pathways such as metabolism of xenobiotics by cytochrome P450 (ko00982) and drug metabolism-other enzymes (ko00983), and some apoptosis related pathways such as p53 signaling pathway (ko04115) and apoptosis pathway (ko04210) were also significantly enriched. These results indicated that the corresponding miRNAs are involved in various processes of shrimp defense

Acknowledgements This research was supported by the National Natural Science Foundation of China (31600321), Guangdong Provincial Natural 159

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

Fig. 7. The top 20 KEGG pathway analysis of novel miRNAs with differential expression. Note: The horizontal coordinates is rich factor. The bigger rich factor the higher enrichment.

Science Foundation (2015A030310438), Program for Scientific Research Start-Up Funds of Guangdong Ocean University (201710566003), Special Program for Outstanding Young Teachers of Guangdong Ocean University (HDYQ2015003). Guangdong Ocean University Student's Platform for Innovation and Entrepreneurship Training.

(2017) 1–8. [5] A. Tassanakajon, K. Somboonwiwat, P. Supungul, S. Tang, Discovery of immune molecules and their crucial functions in shrimp immunity, Fish Shellfish Immunol. 34 (2013) 954–967. [6] V. Matozzo, L. Ballarin, D.M. Pampanin, M.G. Marin, Effects of copper and cadmium exposure on functional responses of hemocytes in the Clam, Tapes philippinarum, Arch Environ Con Tox 41 (2001) 163–170. [7] S.T. Yeh, C.H. Liu, J.C. Chen, Effect of copper sulfate on the immune response and susceptibility to Vibrio alginolyticus in the white shrimp Litopenaeus vannamei, Fish Shellfish Immunol. 17 (2004) 437–446. [8] H. Guo, K. Li, W. Wang, C. Wang, Y. Shen, Effects of copper on hemocyte apoptosis, ROS production, and gene expression in white shrimp Litopenaeus vannamei, Biol. Trace Elem. Res. 179 (2017) 318–326. [9] M.S. Ebert, P.A. Sharp, Roles for microRNAs in conferring robustness to biological processes, Cell 149 (2012) 515–524. [10] D.P. Bartel, MicroRNAs: genomics, biogenesis, mechanism, and function, Cell 116 (2004) 281–297. [11] I.M. Pedersen, G. Cheng, S. Wieland, S. Volinia, C.M. Croce, F.V. Chisari, M. David, Interferon modulation of cellular microRNAs as an antiviral mechanism, Nature 449 (2007) 919–922. [12] V. Ambros, The functions of animal microRNAs, Nature 431 (2004) 350–355. [13] L. Xu, Y. Wang, Y. Xu, L. Wang, L. Zhai, X. Zhu, Y. Gong, S. Ye, L. Liu, Identification and characterization of novel and conserved microRNAs in radish (Raphanus sativus L.) using high-throughput sequencing, Plant Sci 201–202 (2013) 108–114. [14] C. Guo, H. Cui, S. Ni, Y. Yan, Q. Qin, Comprehensive identification and profiling of host miRNAs in response to Singapore grouper iridovirus (SGIV) infection in grouper (Epinephelus coioides), Dev. Comp. Immunol. 52 (2015) 226–235. [15] L. He, G.J. Hannon, MicroRNAs: small RNAs with a big role in gene regulation, Nat. Rev. Genet. 5 (2004) 522–531.

Appendix A. Supplementary data Supplementary data related to this article can be found at http://dx. doi.org/10.1016/j.fsi.2017.12.053. References [1] S.L. Wang, X.R. Xu, Y.X. Sun, J.L. Liu, H.B. Li, Heavy metal pollution in coastal areas of South China: a review, Mar. Pollut. Bull. 76 (2013) 7–15. [2] C.M. Liao, C.F. Chang, C.H. Yeh, S.C. Chen, K.C. Chiang, C.P. Chio, B.Y.H. Chou, L.J. Jou, G.W. Lien, C.M. Lin, H.H. Shen, G.D. Wu, Metal stresses affect the population dynamics of disease transmission in aquaculture species, Aquaculture 257 (2006) 321–332. [3] J.C. Chen, C.H. Lin, Toxicity of copper sulfate for survival, growth, molting and feeding of juveniles of the tiger shrimp, Penaeus monodon, Aquaculture 192 (2001) 55–65. [4] K. Koiwai, R.R.R. Alenton, R. Shiomi, R. Nozaki, H. Kondo, I. Hirono, Two hemocyte sub-populations of kuruma shrimp Marsupenaeus japonicus, Mol. Immunol. 85

160

Fish and Shellfish Immunology 74 (2018) 152–161

H. Guo et al.

[16] L. Ruan, X. Bian, Y. Ji, M. Li, F. Li, X. Yan, Isolation and identification of novel microRNAs from Marsupenaeus japonicus, Fish Shellfish Immunol. 31 (2011) 334–340. [17] T.A. Fehniger, T. Wylie, E. Germino, J.W. Leong, V.J. Magrini, S. Koul, C.R. Keppel, S.E. Schneider, D.C. Koboldt, R.P. Sullivan, M.E. Heinz, S.D. Crosby, R. Nagarajan, G. Ramsingh, D.C. Link, T.J. Ley, E.R. Mardis, Next-generation sequencing identifies the natural killer cell microRNA transcriptome, Genome Res. 20 (2010) 1590–1604. [18] E.R. Mardis, The impact of next-generation sequencing technology on genetics, Trends Genet. 24 (2008) 133–141. [19] E.R. Mardis, Next-generation DNA sequencing methods, Annu. Rev. Genom. Hum. Genet. 9 (2008) 387–402. [20] Q.Y. Xi, Y.Y. Xiong, Y.M. Wang, X. Cheng, Q.E. Qi, G. Shu, S.B. Wang, L.N. Wang, P. Gao, X.T. Zhu, Q.Y. Jiang, Y.L. Zhang, L. Liu, Genome-wide discovery of novel and conserved microRNAs in white shrimp (Litopenaeus vannamei), Mol. Biol. Rep. 42 (2015) 61–69. [21] A. Kozomara, S. Griffiths-Jones, miRBase: integrating microRNA annotation and deep-sequencing data, Nucleic Acids Res. (2010) D152–D157. [22] F. Zhu, Z. Wang, B.Z. Sun, Differential expression of microRNAs in shrimp Marsupenaeus japonicus in response to Vibrio alginolyticus infection, Dev. Comp. Immunol. 55 (2016) 76–79. [23] D.G. Zeng, X.L. Chen, D.X. Xie, Y.Z. Zhao, Q. Yang, H. Wang, Y.M. Li, X.H. Chen, Identification of highly expressed host microRNAs that respond to white spot syndrome virus infection in the Pacific white shrimp Litopenaeus vannamei (Penaeidae), Genet. Mol. Res. 14 (2015) 4818–4828. [24] J. Ou, Q. Meng, Y. Li, Y. Xiu, J. Du, W. Gu, T. Wu, W. Li, Z. Ding, W. Wang, Identification and comparative analysis of the Eriocheir sinensis microRNA transcriptome response to Spiroplasma eriocheiris infection using a deep sequencing approach, Fish Shellfish Immunol. 32 (2012) 345–352. [25] N. Kaewkascholkul, K. Somboonviwat, S. Asakawa, I. Hirono, A. Tassanakajon, K. Somboonwiwat, Shrimp miRNAs regulate innate immune response against white spot syndrome virus infection, Dev. Comp. Immunol. 60 (2016) 191–201. [26] L. Robertson, W. Bray, J. Leung-Trujillo, A. Lawrence, Practical molt staging of Penaeus setiferus and Penaeus stylirostris, J. World Aquacult. Soc. 18 (1987) 180–185. [27] M.R. Friedländer, S.D. Mackowiak, N. Li, W. Chen, N. Rajewsky, miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades, Nucleic Acids Res. 40 (2012) 37–52. [28] H. Guo, C.X. Ye, A.L. Wang, J.A. Xian, S.A. Liao, Y.T. Miao, S.P. Zhang, Trascriptome analysis of the Pacific white shrimp Litopenaeus vannamei exposed to nitrite by RNA-seq, Fish Shellfish Immunol. 35 (2013) 2008–2016. [29] L. Zhou, J. Chen, Z. Li, X. Li, X. Hu, Y. Huang, X. Zhao, C. Liang, Y. Wang, L. Sun, Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma, PLos One 5 (2010) e15224. [30] S. Audic, J.M. Claverie, The significance of digital gene expression profiles, Genome Res. 7 (1997) 986–995. [31] B. Wang, Z. Gan, S. Cai, Z. Wang, D. Yu, Z. Lin, Y. Lu, Z. Wu, J. Jian, Comprehensive identification and profiling of Nile tilapia (Oreochromis niloticus) microRNAs response to Streptococcus agalactiae infection through high-throughput sequencing, Fish Shellfish Immunol. 54 (2016) 93–106. [32] L. Wang, Z. Feng, X. Wang, X. Wang, X. Zhang, DEGseq: an R package for identifying differentially expressed genes from RNA-seq data, Bioinformatics 26 (2010) 136–138. [33] Z. Wang, F. Zhu, MicroRNA-100 is involved in shrimp immune response to white spot syndrome virus (WSSV) and Vibrio alginolyticus infection, Sci. Rep. 7 (2017) 42334. [34] H. Guo, J.A. Xian, A.L. Wang, Analysis of digital gene expression profiling in hemocytes of white shrimp Litopenaeus vannamei under nitrite stress, Fish Shellfish Immunol. 56 (2016) 1–11. [35] X. Chen, Q.B. Li, J. Wang, X. Guo, X.R. Jiang, Z.J. Ren, C.Y. Weng, G. Sun, X.Q. Wang, Y.P. Liu, Identification and characterization of novel amphioxus microRNAs by Solexa sequencing, Genome Biol. 10 (2009) 1–13. [36] C. Cheng, N. Bhardwaj, M. Gerstein, The relationship between the evolution of microRNA targets and the length of their UTRs, BMC Genom. 10 (2009) 431. [37] K.K. Biggar, K.B. Storey, Insight into post-transcriptional gene regulation: stressresponsive microRNAs and their role in the environmental stress survival of tolerant animals, J. Exp. Biol. 218 (2015) 1281–1289. [38] K. Masaomi, K.M. Abul, C. Chao, An intestinal microRNA modulates the homeostatic adaptation to chronic oxidative stress in C. elegans, Aging 8 (2016) 1979–1996. [39] A.P. Gajigan, C. Conaco, A microRNA regulates the response of corals to thermal stress, Mol. Ecol. 26 (2017) 3472–3483. [40] K.K. Biggar, S.F. Kornfeld, Y. Maistrovski, K.B. Storey, MicroRNA regulation in

[41]

[42]

[43] [44]

[45]

[46] [47]

[48] [49] [50]

[51]

[52]

[53]

[54]

[55]

[56]

[57]

[58]

[59]

[60]

[61]

[62]

161

extreme environments: differential expression of microRNAs in the intertidal Snail Littorina littorea during extended periods of freezing and anoxia, Genomics Proteomics Bioinformatics 10 (2012) 302–309. M. Chen, K.B. Storey, Large-scale identification and comparative analysis of miRNA expression profile in the respiratory tree of the sea cucumber Apostichopus japonicus during aestivation, Mar Genom 13 (2014) 39–44. X.P. Li, X.H. Meng, K. Luo, S. Luan, X. Shi, B.X. Cao, J. Kong, The identification of microRNAs involved in the response of Chinese shrimp Fenneropenaeus chinensis to white spot syndrome virus infection, Fish Shellfish Immunol. 68 (2017) 220–231. T.Z. Huang, D.D. Xu, X.B. Zhang, Characterization of host microRNAs that respond to DNA virus infection in a crustacean, BMC Genom. 13 (2012) 1–10. M.B. Ruan, Y.T. Zhao, Z.H. Meng, X.J. Wang, W.C. Yang, Conserved miRNA analysis in Gossypium hirsutum through small RNA sequencing, Genomics 94 (2009) 263–268. T.T. Tan, M. Chen, J.A. Harikrishna, N. Khairuddin, M.I. Mohd Shamsudin, G. Zhang, S. Bhassu, Deep parallel sequencing reveals conserved and novel miRNAs in gill and hepatopancreas of giant freshwater prawn, Fish Shellfish Immunol. 35 (2013) 1061–1069. C.W. Wu, K.B. Storey, Regulation of Smad mediated microRNA transcriptional response in ground squirrels during hibernation, Mol. Cell. Biochem. (2017) 1–11. S. Vidyanand, S. Marepally, S.A. Elliott, S. Baid, V. Lakshmanan, N. Nayyar, D. Bansal, A. Sanchez Alvarado, P.K. Vemula, D. Palakodeti, The miR-124 family of microRNAs is critical for regeneration of the brain and visual system in the planarian Schmidtea mediterranea, Development (Cambridge, England) 144 (2017) 3211–3223. S. Sharma, Immunomodulation: a definitive role of microRNA-142, Dev. Comp. Immunol. 77 (2017) 150–156. J. Luo, P. Chen, W. Xie, F. Wu, MicroRNA-138 inhibits cell proliferation in hepatocellular carcinoma by targeting Sirt1, Oncol. Rep. 38 (2017) 1067–1074. K.J. Peterson, M.R. Dietrich, M.A. McPeek, MicroRNAs and metazoan macroevolution: insights into canalization, complexity, and the Cambrian explosion, Bioessays 31 (2009) 736–747. J. Wen, B.J. Parker, G.F. Weiller, In Silico identification and characterization of mRNA-like noncoding transcripts in Medicago truncatula, Silico Biol. 7 (2007) 485–505. X. Sun, Q.-h. Liu, B. Yang, J. Huang, Differential expression of microRNAs of Litopenaeus vannamei in response to different virulence WSSV infection, Fish Shellfish Immunol. 58 (2016) 18–23. F.F. Fagutao, T. Koyama, A. Kaizu, T. Saito-Taki, H. Kondo, T. Aoki, I. Hirono, Increased bacterial load in shrimp hemolymph in the absence of prophenoloxidase, FEBS J. 276 (2009) 5298–5306. S. Radha, P. Mullainadhan, M. Arumugam, Detection of two distinct types of hemolymphatic prophenoloxidase and their differential responses in the black tiger shrimp, Penaeus monodon, upon infection by white spot syndrome virus, Aquaculture 376 (2013) 76–84. S. Boonrawd, R. Mani, S. Ponprateep, P. Supungul, P. Masrinoul, A. Tassanakajon, V. Rimphanitchayakit, Characterization of PmSpӓtzle 1 from the black tiger shrimp Peneaus monodon, Fish Shellfish Immunol. 65 (2017) 88–95. Y.R. Shi, M. Jin, F.T. Ma, Y. Huang, X. Huang, J.L. Feng, L.L. Zhao, Y.H. Chen, Q. Ren, Involvement of Relish gene from Macrobrachium rosenbergii in the expression of anti-microbial peptides, Dev. Comp. Immunol. 52 (2015) 236–244. Y. Ren, H. Zhang, B. Pan, C. Yan, A Kazal-type serine proteinase inhibitor from Cyclina sinensis is involved in immune response and signal pathway initiation, Fish Shellfish Immunol. 47 (2015) 110–116. Y. Qiao, Y. Mao, J. Wang, R.N. Chen, L.B. Zheng, Y.Q. Su, J. Chen, W.Q. Zheng, Analysis of liver and gill miRNAs of Larimichthys crocea against Cryptocryon irritans challenge, Fish Shellfish Immunol. 59 (2016) 484–491. H. Chen, L.S. Xin, X.R. Song, L. Wang, W.L. Wang, Z.Q. Liu, H. Zhang, L.L. Wang, Z. Zhou, L.M. Qiu, L.S. Song, A norepinephrine-responsive miRNA directly promotes CgHSP90AA1 expression in oyster haemocytes during desiccation, Fish Shellfish Immunol. 64 (2017) 297–307. H. Wang, J.J. Kim, E. Denou, A. Gallagher, D.J. Thornton, M.S. Shajib, L. Xia, J.D. Schertzer, R.K. Grencis, D.J. Philpott, New role of Nod proteins in regulation of intestinal goblet cell response in the context of innate host defense in an enteric parasite infection, Infect. Immun. 84 (2016) 275–285. A. Tawiah, F. Moreau, K. Chadee, High MUC2 production in goblet cells induces ER stress and exhibit increase susceptibility to apoptosis, Faseb. J. 28 (2014) 9601–9604. Z. Jiang, T.J. Applegate, A.C. Lossie, Cloning, annotation and developmental expression of the chicken intestinal MUC2 gene, PLos One 8 (2013) e53781.