Dynamic mRNA and miRNA expression analysis in response to hypoxia and reoxygenation in the blunt snout bream (Megalobrama amblycephala)

Adaptation to hypoxia is a complex process involving various pathways and regulation mechanisms. A better understanding of the genetic influence on these mechanisms could permit selection for hypoxia-sensitive fish. To aid this understanding, an integrated analysis of miRNA and mRNA expression was performed in Megalobrama amblycephala under four acute hypoxia and reoxygenation stages. A number of significantly differentially-expressed miRNAs and genes associated with oxidative stress were identified, and their functional characteristics were revealed by GO function and KEGG pathway analysis. They were found to be involved in HIF-1 pathways known to affect energy metabolism and apoptosis. MiRNA-mRNA interaction pairs were detected from comparison of expression between the four different stages. The function annotation results also showed that many miRNA-mRNA interaction pairs were likely to be involved in regulating hypoxia stress. As a unique resource for gene expression and regulation during hypoxia and reoxygenation, this study could provide a starting point for further studies to better understand the genetic background of hypoxia stress.

Cyprinus carpio 17 . Its sensitivity to hypoxia is such that even a short period (<2 h) of hypoxia (<0.5 mg/L) at room temperature can be lethal 18 . Although this sensitivity to hypoxia has implications for the welfare and culture of M. amblycephala, it could allow this species to be used as a model organism to explore the molecular mechanisms of acute hypoxia and reoxygenation via the transcriptome.
Next-generation sequencing (NGS)-based RNA sequencing (RNA-seq) methods of transcriptome analysis not only allow sequences to be acquired for gene discovery but also help to identify which transcripts are involved in hypoxic stress [19][20][21] . Recent advances in RNA-seq have generated an unprecedented global view of the transcriptome and have also provided a more efficient method to explore the transcriptional landscape 22,23 . MicroRNA (miRNA) is one of the epigenetic mechanisms which can regulate gene expression in a tissue-specific manner 24 . MiRNAs have been identified in a wide variety of organisms and are currently being investigated in laboratory fish models such as the zebrafish Danio rerio 25 , the Japanese medaka Oryzias latipes 26 and marine medaka Oryzias melastigma 27 . Integrated analysis of mRNA and miRNA expression profiles has been performed in the zebrafish for optic nerve regeneration, in M. amblycephala in response to intermuscular bone development and in darkbarbel catfish (Pelteobagrus vachelli) in response to hypoxia [28][29][30] . A possible next step for research could therefore be to investigate whether hypoxia-induced dysregulation of miRNA causes oxidative tress impairment, and if so to identify the mechanism. A potentially more reliable method than in vivo trials for predicting miRNA-mRNA interactions is to integrate real mRNA and miRNA transcriptomic data into in silico target predictions 31 .
We therefore used NGS to determine mRNA-seq and miRNA-seq in the livers of both control and hypoxia-treated M. amblycephala, to investigate the molecular mechanisms of hypoxia adaptation. We identified hypoxia time specific expressed mRNAs and miRNAs to further understand the molecular mechanisms of hypoxia adaptation. The study also found various interaction networks and regulatory modes of mRNAs and miRNAs, based on the integrated analysis of miRNA and mRNA expression profiles. Finding the molecular mechanisms of hypoxia adaptation in fish will not only help us to understand the evolution of the hypoxia-signalling pathway but will also inform breeding of hypoxia-tolerant fish strains.

Results and Discussion
Assembly and annotation of reference transcriptome. To obtain a reference transcriptome for M. amblycephala in response to hypoxia and reoxygenation, a RNA-seq library was constructed using RNA from samples of four stages. A total of 425,901,164 raw reads were generated using high-throughput sequencing. After quality control, approximately 416,787,182 high-quality reads with a Q20 percentage of 96.51% and GC percentage of 46.08% were available for analysis (Table 1). A total of 293,928 contigs were generated using the Trinity de novo assembler, with an average length of 566 bp and an N50 of 667 bp. A total of 279,747 unigenes were generated with an average length of 534 bp ( Table 1). The length distribution of these contigs is shown in Supplementary Figure S1. A total of 89,118 unigenes were successfully annotated via alignment to reference databases. The NR, Swiss-Prot, KEGG and KOG databases were used to annotate 70,747 (79.39%), 58,498 (65.64%), 55,606 (62.40%) and 40,166 (45.07%) unigenes, respectively (Table 1). Annotation was based on the NCBI NR database, E-value distribution, similarity distribution and species distribution of the result of NR annotation (Supplementary Figure S2). All the data are available at the NCBI SRA database (SRP100308).
To identify the functional genes of M. amblycephala in response to hypoxia and reoxygenation, 12 gene expression profiling libraries were generated from the four stages, with three biological replicates. The mapping results for the 12 libraries are shown in the Table 2, and the distributions of unique reads were compared to indicate gene coverage. Similar coverage was obtained for all 12 libraries. Principal component analysis showed that the biological replicates had very similar expression levels, suggesting good reproducibility of the method (Supplementary Figure S3). A total of 7,756 differential expressed genes (DEGs) were found between hypoxia 3 h, hypoxia 24 h, reoxygenation 3 h and normoxia conditions. Of these DEGs, 4658, 2500 and 2635 genes were expressed in response to hypoxia 3 h, hypoxia 24 h and reoxygenation 3 h, respectively (Fig. 1A). The Venn diagram (Fig. 1B) of the DEGs shows that most genes were not expressed in all hypoxia stages. This suggests that the mechanisms and pathways involved in response to hypoxia are different during hypoxia and post-hypoxia recovery (Fig. 1B). Hierarchical clustering of partial differentially-expressed mRNAs was seen between different stages (Fig. 1C). Annotation revealed hundreds of hypoxia-related genes in the transcriptome (Supplementary  Functional analysis of DEGs during hypoxia and reoxygenation. The liver in fish has relatively high metabolic oxygen demands 32 . Previous studies of expression profiling in fish have found that a cluster of genes involved in oxygen utilization are differentially expressed in response to hypoxia [33][34][35][36] . An increased hemoglobin mRNA expression levels has been found during hypoxia acclimation in liver of M. amblycephala, which was similar with previous study in muscle of this specie 37 . It was beneficial for utilization of circulating oxygen during hypoxia as that could facilitate the diffusion of oxygen. The electron transport chain in mitochondria is a major site of ROS production. Hypoxic stress partially inhibits mitochondrial electron transport, leading to redox changes in the electron carriers and increasing ROS production at complex III 38,39 . A compromised antioxidant defence or inhibition of electron flow can disrupt the delicate balance between antioxidant defence and ROS production 40 , thus leaving M. amblycephala susceptible to ROS damage during hypoxia. Interestingly, superoxide dismutase (SOD), catalase (CAT) and nuclear factor erythroid 2-related factor 2 (Nrf2) was significantly higher in M. amblycephala in response to reoxygenation 3 h, previous study have confrimed that the Nrf2 in modulating the activation of antioxidant defences at both transcriptional and functional levels under reoxygenation in Rhamdia quelen 41 , Cyprinus carpio 42 and Pelteobagrus vachelli 43 . Although the antioxidant enzyme levels could only reflect the antioxidant capacity rather than the status of oxidative stress, previous study confrimed that ROS still significant increased in M. amblycephala during the recovery from severe acute hypoxia (DO: 1.0 ± 0.2 mg/L) treatment 44 , which supports the hypothesis that anticipatory preparation takes place in order to deal with the encountered oxidative stress during the recovery from hypoxia as proposed by M. Hermes-Lima.
Enrichment analysis is an effective way to identify those KEGG pathways which frequently occur in a tissue, against a background of the whole body transcriptome 45,46 . A total of 532 DEGs were mapped into 108 KEGG pathways in response to hypoxia and reoxygenation (Supplementary Table S2). By performing the KEGG pathway analyses, top 20 significantly changed pathways were identified in response to hypoxia and reoxygenation, respectively ( Figure S4). Among these top 20 pathways, "Signal transduction" were the most commonly represented subclasses in hypoxia 3 h vs. normoxia group. "HIF-1 signaling pathway", "Insulin signaling pathway", "AMPK signaling pathway", "P53 signaling pathway" and "FoxO signaling pathway" are included in the "Signal transduction" subclass, suggesting that a number of other transcription factors are activated either directly or indirectly in response to acute hypoxia. "Immune system" were the most commonly represented subclasses in hypoxia 24 h vs. normoxia group. "Viral myocarditis", "Viral carcinogenesis", "Phagosome", "Complement and coagulation cascades" are included in the "Immune system" subclass. The present results imply that our prolonged hypoxia treatment affected immunity-related genes involved in these pathways in M. amblycephala. "Glycolysis/ Gluconeogenesis", "Fructose and mannose metabolism", "Galactose metabolism", "Pentose phosphate pathway", "Starch and sucrose metabolism", "Fatty acid metabolism" and "Biosynthesis of amino acids" have also been significantly enriched in reoxygenation 3 h vs. normoxia group. The enrichment of "Substance metabolism" subclasses may be necessary to cope with the energy demand balance from the liver of M. amblycephala in respond to reoxygenation 3 h.
While the hypoxia-inducible factor (HIF-1) plays a major role in controlling the ubiquitous transcriptional response to hypoxia, we focused our analyses on genes likely to be relevant to hypoxia 3 h and found that the HIF-1 signalling pathway co-operated with the AMP-activated protein kinase (AMPK) pathway to regulate and control hypoxia stress (Supplementary Figure S5). HIF-1α is rapidly broken down by prolyl hydroxylases in normoxia, but is stabilised in hypoxia [47][48][49][50] . This suggests a model of hypoxia-induced regulation of the HIF-1α transcript in fish, which could be similar to typical HIF-1α transcriptional regulation in mammals 51 53 . AMPK is activated by increases in the cellular AMP:ATP ratio caused by hypoxia stresses that either reduce ATP production or ATP-producing catabolic pathways such as fatty acid oxidation and glycolysis 54,55 , suggesting that anaerobic metabolism has probably taken the place of aerobic metabolism during hypoxia 3 h since M. amblycephala suffocation point ranges from 0.64 to 0.35 mg/L 56 . Aerobic metabolism is less efficient than aerobic respiration when producing energy, and a considerable increase in glycolytic flux is needed to avoid a detrimental decrease of cellular ATP [57][58][59] . The KEGG analysis revealed that some pathways related to immune response were significantly down-regulated in response hypoxia 24 h, it was suggested that hypoxia stress can modulate and compromise the innate immune response as previous shown in M. amblycephala 60 . It is noteworthy that "Oxidative phosphorylation" was found to be the most frequently represented subclasses in M. amblycephala in response to reoxygenation 3 h, suggesting increase total cellular ATP content by induced glycolysis and mitochondria oxidative phosphorylation in fish for reoxygenation. In addition, similar with previous study in M. amblycephala 44 , some important immunity-related pathway have also been significantly enriched, suggesting that the difference of energy metabolism pattern and immune response is the main biological process of hypoxia and recover.
Sequence and expression profiling of hypoxia-related microRNAs.   Table S4). Using a |log 2 FC| of ≥1 and reads ≥100 as the cut-off, we identified the following differentially-expressed miRNAs in adjacent pairwise comparisons ( Fig. 2A): in normoxia-vs-hypoxia 3 h there were 132 (40 up-regulated and 92 down-regulated); in normoxia-vs-hypoxia 24 h there were 120 (69 up-regulated and 51 down-regulated); and in normoxia-vs-reoxygnation 3 h there were 174 (108 up-regulated and 66 down-regulated) miRNAs. Through Venn diagram analysis, 33 overlapping differentially-expressed miRNAs were identified in the three pairwise and three nonadjacent pairwise comparisons (Fig. 2B). Hierarchical clustering of partial differentially-expressed miRNAs was seen between different stages (Fig. 2C).
In this study, we found several miRNAs that were significantly differentially expressed in response to hypoxia. For example, miR-9 and miR-124 are both known to play a crucial role in neurogenesis and neuronal development, particularly in the regulation of neural differentiation, proliferation and cell migration [61][62][63] . Both miR-143 and miR-101 can directly target the key glycolytic enzyme hexokinase under hypoxic conditions in tissue cells 64,65 , as we known, hexokinase is the first key regulatory enzyme of the glycolytic pathway 66 , previous study show that hexokinase could accelerates the rate of glycolysis to produce ATP anaerobically in Gadus morhua in response to hypoxia stress 36 , which was consitent with our results of this study in M. amblycephala (Fig. 3). In addition, miR-NAs implicated in the regulation of apoptosis include miR-210, miR-21, miR-17 and miR-29, which were found to be expressed in the current study. Of these, miR-210 and miR-21 exert broad pleiotropic effects by targeting genes involved in mitochondrial function, apoptosis, cell cycle arrest and cell survival 67,68 , while miR-29 is thought to be pro-apoptotic, suppressing tumour growth in mice 69,70 . Those miRNA above were predicted target gene including p53, caspase-3 and caspase-9, were significant increased in liver of M. amblycephala in response to hypoxia (Fig. 3), The tumor suppressor protein p53 is a universal sensor of environmental stress, it plays a critical role in promoting the survival or death of cells exposed to agents that cause DNA damage 71 , caspase 3 and caspase-9 play a pivotal role in mitochondria apoptosis pathway 72 , suggesting apoptosis pathway in M. amblycephala as a hypoxic sensitive fish was induced by hypoxia stress. Therefore, some of the miRNAs that are differentially expressed between hypoxia and normoxia stages are likely to play important roles in the physiological response. Quantitative analysis of miRNA and mRNA expression. To validate the RNA-seq results, 18 genes that showed a high level expression (or are known to be important in the stress response) were selected for qRT-PCR analysis. Same expression trends were found between qRT-PCR and the Illumina data (Fig. 3). Eight mRNAs were manually selected as representatives for their potential roles in hypoxia response according to their annotations and their potential relationship with four hypoxia-responsive miRNAs. These genes encode AMP-activated protein kinase (miR-92a-2-5p-AMPK), lactate dehydrogenase (miR-92a-2-5p-LDH), NADH dehydrogenase subunit 1 (miR-92a-2-5p-ND1), acetyl-coenzyme A synthetase (miR-222b-AcCoA), lactate dehydrogenase (miR-17-3p-LDH), pyruvate kinase (miR-17-3p-PK), 1-phosphatidylinositol 4,5-bisphosphate phosphodiesterase (miR-17-3p-PLCG), cytochrome coxidase subunit (miR-132b-CCO). A similar result was seen in the validation of the miRNA-mRNA interaction pairs: four interaction pairs had the same expression pattern as the sequencing data (Fig. 4). We believe the results confirm the credibility of the molecular resources and sequencing data identified in our study.
The complexity of the miRNA-mRNA interaction pairs involved in specific biological processes has been an exciting and challenging field of study 73 , Due to the limited genetic background of M. amblycephala and that the present ESTs may be the noncoding RNA, our study describes here the interaction between 65 genes with annotation and 8 miRNAs, those miRNA-mRNA interaction network during early hypoxia 3 h help us understand the roles for specific miRNAs or miRNA-mRNA interactions (Fig. 5).

Conclusion
This study is the first integrated analysis of miRNA and mRNA expression during hypoxia and reoxygenation in a fish species. Putative miRNAs and genes associated with hypoxia were identified in the study. We found the obvious negative correlation between miRNA and their target mRNAs, providing several miRNA-mRNA interaction networks in M. amblycephala in response to hypoxia. This study generated fundamental molecular resources across four hypoxia and reoxygenation stages, which can be used to further understand the molecular processes of hypoxia stress. Fouther analysis of miRNA-mRNA regulatory network from other tissues or cell types by bioinformatics analysis will be required.

Materials and Methods
Ethics statement. All animals and experiments were conducted in accordance with the "Guidelines for Experimental Animals" of the Ministry of Science and Technology (Beijing, China). All efforts were made to minimize suffering. All experimental procedures involving fish were approved by the institution animal care and use committee of the Chinese Academy of Fishery Sciences.
Sample collection. Two hundred of healthy M. amblycephala juveniles with a mean weight of 50.56 ± 5.15 g were obtained from the Freshwater Fisheries Research Center, Chinese Academy of Fishery Sciences, China. The fish were immediately transferred to the aquatic laboratory and kept in six 500-L fiberglass tanks (n = 30 fish/ tank). During acclimation, each tank was filled with pre-aerated municipal water kept under the following conditions: 24.6 ± 0.5 °C, pH 8.2 ± 0.08, dissolved oxygen 6.5 ± 0.2 mg/L, total ammonia-nitrogen 0.08-0.09 mg/L, with natural photoperiod. The control group was maintained under normoxic conditions (6.5 ± 0.2 mg/L). In the treatment tanks, hypoxic conditions (1.5 ± 0.1 mg/L dissolved oxygen) were maintained by nitrogen manipulation 74 . Fish were randomly allocated to one of four treatment groups: hypoxia for 0 h (normoxia, control group), hypoxia for 3 h, hypoxia for 24 h and hypoxia for 24 h + 3 h recovery in normoxia. Liver samples were obtained RNA isolation, cDNA library construction and sequencing. Total RNA was extracted from the tissue using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's directions. RNA samples were digested using DNase I to remove potential genomic DNA. The integrity and size distribution were assessed using a Bioanalyzer 2100 with RNA 6000 Nano Labchips (Agilent technologies, Santa Clara, CA, USA), and all samples were standardised to 500 ng/μl. Three sequencing libraries and one miRNA transcriptome library were constructed for each hypoxia and reoxygenation condition. Next, 12 sequencing libraries were constructed using a TruSeq ™ RNA Sample Preparation Kit (Illumina) according to the manufacturer's directions. The next step was KAPA quantification and dilution, after which the library was sequenced on an Illumina HiSeq. 2500 with 125-bp paired-end reads. After removing adaptor sequences, ambiguous N nucleotides (those an N ratio <5%) and low quality sequences (reads with <50% bases of quality value), the remaining clean reads were assembled using Trinity software 75 . This was performed as described for a de novo transcriptome assembly without a reference genome. The process generated the reference sequences, including several EST sequences available in the public domain for comparative transcriptome study. All data were made available in the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra/?term=SRP100308).
Functional annotation of assembled contigs. Functional annotation of the assembled reference sequences was performed by homology searches against the NCBI non-redundant (NR) protein database 76 , The Universal Protein Resource (UniProt-SwissProt) database 77 and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database 78,79 . The assembled transcriptome contigs were subjected to similarity searches against the NCBI NR protein database using BLASTx with an e-value cut-off of 1e-10. A gene name and description was assigned to each contig based on the BLASTx hit with the highest score.
Differential gene expression analysis. The cleaned reads were mapped to the assembled reference transcriptome using Bowtie 80 , and about 70% of the reads for each sample could be mapped to the reference (Table 1). Gene and isoform abundances were quantified using RSEM according to the Trinity-assembled transcriptome. Gene expression was measured as fragments per kilobase of exon per million reads mapped (FPKM). Finally, we normalised the expression levels in each sample using edgeR and obtained the differentially-expressed transcripts using pairwise comparisons 81 . The following threshold values were used to judge the significance of DEGs: a |log 2 (fold change)| of ≥1 and a false discovery rate of ≤0.05. GO analysis was conducted on the assembled transcriptome using InterProScan (http://www.ebi.ac.uk/Tools/pfa/iprscan/) and integrated protein databases with default parameters. The GO terms associated with transcriptome contigs were then obtained for describing their biological processes, molecular functions and cellular components 32 . For pathway enrichment analysis, all DEGs were mapped to terms in the KEGG database and searched for significantly enriched KEGG terms compared to the whole transcriptome background 32 . GO and KEGG were performed using the ultrageometric test to identify which DEGs were significantly enriched in GO terms (P-value ≤ 0.05) and KEGG pathways (q-value ≤ 0.05) compared with the whole transcriptome background. MiRNA sequencing and differential expression analysis. Four small RNA libraries were constructed and sequenced as previously described 82 . The Rfam database (ftp://ftp.sanger.ac.uk/pub/databases/Rfam/9.1/) database and the GenBank noncoding RNA database (http://blast.ncbi.nlm.nih.gov/) were searched using BLAST for the high-quality reads obtained. Once the rRNA, tRNA, snRNA and other ncRNA sequences were annotated, they were then aligned to mRNA exons and introns to screen and remove degraded fragments. We also mapped selected sequences to the reference transcriptome (with a tolerance of one mismatch in the seed sequence) to analyse their expression and distribution on the genome. To determine novelty and species specificity, the remaining clean reads were aligned against known deuterostome pre-miRNAs in miRbase 20.0.
Differences in miRNA expression between the hypoxia and reoxygenation stages were determined as |log 2 (fold change)| >1, and significance was set at a P-value of ≤0.05. Subsequently, differences in gene expression which were complementary with corresponding miRNA expression values were selected and analysed with both Targetscan (http://www.targetscan.org/) and MiRanda (http://www.microrna.org) to predict the miRNA target 83 . Finally, miRNA-mRNA pairs were detected and their interaction analysed based on the target prediction, function annotation and negative regulation mechanism of mRNA and miRNA. The resulting miRNA-mRNA interaction networks were displayed using visualising maps.
Quantitative PCR for miRNA and mRNA expression. The sequencing data were validated as follows: 18 mRNAs which showed significant increases or decreases in expression among "hypoxia 3 h vs normoxia", "hypoxia 24 h vs normoxia" and "reoxygenation 3 h vs normoxia" stage were identified, along with four miRNA-mRNA interaction pairs from the paired comparison of normoxia and hypoxia for 3 h. Specific qRT-PCR primers and microRNA stem-loop RT were used in this study (Supplementary Table S6). High quality total RNA, miRNAs and mRNAs were obtained and then reverse transcribed using a PrimeScript RT reagent kit with gDNA Eraser (TaKaRa, Japan). Quantitative real-time PCR (RT qPCR) analyses of the miRNAs and mRNAs were performed using SYBR Green PCR Master Mix (TaKaRa) according to the manufacturer's directions. The internal RT qPCR controls for the miRNA and mRNA were M. amblycephala U6 snRNA and β-actin, respectively. All PCRs were performed in triplicate. Relative expression levels were measured in terms of threshold cycle value (Ct) and were normalised using the equation 2 −ΔΔCt 84 .