Identification of microRNAs in Silver Carp (Hypophthalmichthys molitrix) Response to Hypoxia Stress

Simple Summary Hypoxia stress is one of the main problems in silver carp (Hypophthalmichthys molitrix) culture. Severe hypoxia stress can lead to damage and even death of silver carp. Therefore, it is very important to explore how silver carp adapt to and respond to hypoxia stress. MicroRNAs play an important role in a series of important life activities in organisms. In this study, the differentiallyexpressed miRNAs were screened from a mixed pool of liver, brain, heart and gill of silver carp under different levels of hypoxia stress by high-throughput sequencing. Our findings provided new insights to further study the miRNA regulatory mechanism and molecular characteristics of anoxic response in silver carp. Abstract Hypoxia is one of the serious stresses in fish culture, which can lead to physical and morphological changes, and cause injury and even death to fish. Silver carp (Hypophthalmichthys molitrix) is an important economic fish and widely distributed in China. MicroRNA is a kind of endogenous non-coding single-stranded small RNA, which is involved in cell development, and immune response and gene expression regulation. In this study, silver carp were kept in the closed containers for hypoxia treatment by spontaneous oxygen consumption. The samples of heart, brain, liver and gill were collected, and the total RNAs extracted separately from the four tissues were mixed in equal amounts according to the concentration. Afterwards, the RNA pool was constructed for high-throughput sequencing, and based on the small RNA sequencing, the differentially expressed microRNAs were identified. Furthermore, their target gene prediction and enrichment analyses were carried out. The results showed that a total of 229 known miRNAs and 391 putative novel miRNAs were identified, which provided valuable resources for further study on the regulatory mechanism of miRNAs in silver carp under hypoxia stress. The authors verified 16 differentially expressed miRNAs by qRT-PCR, and the results were consistent with small RNA sequencing (sRNA-seq). The predicted target genes number of differentially expressed miRNAs was 25,146. GO and KEGG functional enrichment analysis showed that these target genes were mainly involved in the adaption of hypoxia stress in silver carp through biological regulation, catalytic activity and apoptosis. This study provides references for further study of interaction between miRNAs and target genes, and the basic data for the response mechanism under hypoxia stress in silver carp.


Introduction
The phenomenon of dissolved oxygen (DO) under 2 mg/L in the water environment is called hypoxia [1]. Hypoxia is an increasing threat factor to aquaculture and can lead to physical and morphological changes, and cause injury and even death to fish [2,3]. In the process of aquaculture, due to high temperature and high cultivation density, the dissolved oxygen in the water is often lower than the normal level, posing a serious threat to the cultured animals [4,5]. Like other fish, silver carp are prone to suffer from hypoxia stress [6]. Hypoxia can cause apoptosis of hepatocytes and compromise the central nervous system and immune system of fish, leading to death [7]. So far, there have been many studies on the effects of hypoxia on fish species, such as the function of zebrafish (Danio rerio) miRNA in cellular adaptation under hypoxia [8].
Silver carp, Hypophthalmichthys molitrix, is widely distributed in various water bodies in China and is a very popular food because of its fresh meat and rich nutrition. In past years, silver carp rapidly became one of dominant species in aquaculture due to its superior characteristics such as a short food chain, strong disease resistance and low price [9]. Furthermore, it also has a positive effect on improving water quality and plays a vital role in water purification fisheries because it feeds on plankton in water and controls the blue-green algae bloom [10]. However, the hypoxia tolerance of silver carp is very poor, and hypoxia stress usually leads to physical damage and even death due to high-density culture and transportation [11,12]. Hypoxia has become one of the main reasons restricting the survival rate of cultured silver carp. Therefore, studying the hypoxia resistance of silver carp is not only of great significance to the healthy development of silver carp aquaculture, but also has a deeper understanding of the molecular mechanism of hypoxia tolerance of all fish.
The microribonucleic acids (miRNAs) are a kind of endogenous non-coding singlestranded small RNAs with approximately 18-23 nucleotides in length, which are bound to the 3 untranslated region (UTR) of the target gene, thereby repressing and degrading the translation of the target gene mRNA, and playing an important regulatory role [13][14][15]. MiRNAs exist widely in animals, and participate in a series of important life activities of organisms, such as growth and development, apoptosis and immunity [16]. In the largemouth bass (Micropterus salmoides), the miRNA-mRNA was comprehensively analyzed by high-throughput sequencing technology and 13 differential expressions miRNAs of liver involved in glucose and lipid metabolism were identified under acute hypoxia [17]. The miR-462/miR-731 cluster of zebrafish was significantly up-regulated under hypoxia stress, and cell-cycle progression was blocked to inhibit cell proliferation and induce apoptosis [8]. In nile tilapia (Oreochromis niloticus), miR-204, as an endogenous regulator of vascular endothelial growth factor (VEGF) expression, was inhibited under hypoxia stress, and the VEGF expression level was significantly up-regulated [18]. Recently, a total of 324 miRNAs, including 309 conserved miRNAs and 15 novel miRNAs, were identified by deep sequencing of tested tissues (kidney, spleen, muscle and liver) from the blunt snout bream (Megalobrama amblycephala) [19]. High-throughput sequencing technology of miRNAs is an effective method to investigate the expression of miRNAs [20]. Some recent studies on aquatic animals showed that miRNAs were also involved in the process of adapting to changes in the external environment [21]. MiRNA could respond to low salinity stress in the oyster (Crassostrea gigas) [22], heat stress of the sea cucumber (Apostichopus japonicus) [20], and the adaptation of Pacific whiteleg shrimp (Litopenaeus Vannamei) to hypoxia stress [21]. However, the roles of miRNAs under hypoxia stress have not been discussed yet in silver carp.
In this study, we employed a high-throughput sequencing technique to determine the expression of miRNAs in silver carp under hypoxia stress. Gene ontology and KEGG pathway analyses of putative target genes were also carried out. In addition, we also analyzed differentially expressed miRNAs. The data obtained help to clarify the regulatory role of miRNAs in the hypoxic stress response of silver carp, and provide new insights for the study of miRNAs regulation and molecular adaptation mechanisms of silver carp under hypoxia stress.

Ethics Approval
All experimental procedures were conducted according to guidelines of the appropriate Animal Experimental Ethical Inspection of Laboratory Animal Centre, Yangtze River Fisheries Research Institute, Chinese Academy of Fishery Sciences (approval number 2020098).

Experimental Fish and Sample Collection
Silver carp (53.64 ± 4.84 g) in this study were taken from Yaowan Experimental Farm, Yangtze River Fishery Research Institute, Chinese Academy of Fisheries Sciences. The silver carp were temporarily cultured in a tank for a week. Before the experiment, the experimental fish were transferred to the 52.5 L tank for 24 h, and the normal aeration was carried out to maintain the oxygen concentration in the water. A total of 120 healthy silver carp were selected and divided into a normoxia group (T0) and three hypoxia stress groups (T1, T2, T3). They were placed in white water tank (50 cm × 35 cm × 30 cm), and 10 fish were placed in each tank, with 3 replicates in each group. Normal aeration was maintained in the normoxic group (all the fish breathed normally, recorded as T0, dissolved oxygen 6.45 mg/L); hypoxia stress experimental group stopped aeration, sealed water tank with plastic film, respectively. Silver carp gasping for air period (most fish tried to breathe directly through their mouths, recorded as T1, 4 h after the beginning of the experiment, dissolved oxygen 0.76 mg/L), semi-asphyxia period (half of fish lost balance, recorded as T2, 5 h after the beginning of the experiment, dissolved oxygen 0.58 mg/L), asphyxiation period (half of fish sank without the rhythmical opening and closing of the gill flaps, recorded as T3, 6 h after the beginning of the experiment, dissolved oxygen 0.21 mg/L). The fish were sampled, and three silver carp were collected in each water tank at each experiment point after MS-222 (100 mg/L) anesthesia. The heart, liver, brain and gill were kept in 2 mL cryopreservation tube and put into liquid nitrogen, and then stored at −80 • C.

Total RNA Extraction, Library Construction and Small RNA Sequencing
Total RNA was extracted from the heart, liver, brain and gill by Trizol Reagent (Invitrogen). RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer ® spectrophotometer (IMPLEN, München, Germany). RNA concentration was measured using Qubit ® RNA Assay Kit in Qubit ® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Finally, the total RNA extracted from four tissues of different groups was mixed in equal amounts to construct the RNA pools of the normoxic group and the hypoxic stress groups.
The sequencing libraries were further constructed using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® (NEB, San Diego, CA, USA) following manufacturer's instructions. Briefly, RNA adapters were ligated to 3 and 5 end of RNA followed by cDNA synthesis and PCR amplification. The cDNA library was separated by PAGE gel, and small RNAs of 18-40 bp were isolated and purified. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions ( Figure S3). After the clusters were generated, the libraries were sequenced on the Illumina Noveseq platform to generate 50 bp single-end reading codes.

Identification and Analysis of Differentially Expressed miRNAs
After sequencing, clean reads were obtained by removing reads containing adapter, ploy-N and low quality reads (reads having >50%, bases with quality score ≤5) from raw data. Only the remaining reads were regarded as clean reads and used in following analysis. Finally, clean reads were used for subsequent data analysis. Clean reads were compared to the miRNA sequences specified in the miRbase(v22) database by bowtie2 (2.2.2) software to identify known miRNAs. Using DESeq2 to assess differentially expressed miRNAs, miR-NAs with |log2(fold change)| > 1 and adjusted p-value < 0.05 were considered significant ( Figure S1). Target genes of miRNAs were predicted by MiRanda (v3.3a) and qTar software. The final target set was the intersection of the two tools. Subsequently, GO and KEGG enrichment analysis was performed on the target genes of miRNAs to determine the main biochemical metabolic pathways and signal transduction pathways involved. Only GO terms with a Bonferroni-corrected p value ≤ 0.05 and pathways with a FDR ≤ 0.05 were regarded as significantly enriched.

Real-Time PCR Validation
RNA samples from the normal oxygen group and hypoxia stress group were subjected to stem-loop qPCR analysis. MiRNA 1st Strand cDNA Synthesis Kit (Vazyme, Nanjing, China) was used to remove genomic DNA from total RNA and reverse transcribe into cDNA. Quant Studio 5 system was used for amplification reaction, and the amplification system is as follows: 2×ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) 10 µL, forward primer (10 µmol/L) 0.4 µL, reverse primer (10 µmol/L) 0.4 µL, template cDNA 1 µL, ddH 2 O 8.2 µL. The qRT-PCR conditions were carried out as follows: initial denaturation at 95 • C for 30 s, followed by 40 cycles of 30 s denaturation at 95 • C, annealing at 60 • C for 30 s and extension at 72 • C for 20 s. The specificity of PCR amplification was verified by melting curve. U6 was used as an internal control and Premier 5.0 software was used for primer design; stem-loop primer sequence and qPCR primer sequence are listed in Table 1. In addition, 16 differentially expressed miRNAs were randomly selected and their relative expression levels were determined by qRT-PCR to verify the reliability of Illumina Noveseq high-throughput sequencing results. The relative expression levels of miRNAs were calculated by 2 −∆∆Ct method [14]. SPSS software, version 26.0 (IBM Corp, Armonk, NY, USA) was used for the statistical analysis.    Figure S2). The data quality of the small RNA library was listed in Table 2. Note: Sample raw sequencing data yield and quality statistics. Q20 represents a 1% chance that the base will be incorrectly determined; Q30 represents a 1‰ chance that the base will be incorrectly determined.
The clean reads of each sample were screened for small RNA in a certain length range, and it was found that the sequence length was mainly concentrated between 21-23 nt, of which 22 nt accounted for the largest proportion ( Figure 1). A total of 288, 286, 249 and 280 known miRNAs were identified in T0, T1, T2, and T3 group, respectively, while 305, 279, 165 and 298 novel miRNAs were predicted respectively (Table S1). The expression levels of these miRNAs are different, but several miRNAs, such as dre-miR-100-5p, dre-let-7e, dre-miR-101a, dre-miR-143 and dre-miR-146a, were all highly expressed in different levels of hypoxia stress (Table 3). According to the exon, intron and repeat sequence information in genome annotation, and the sRNA information in mirBase and Rfam databases, the small RNA was classified and annotated (Table 4). and 280 known miRNAs were identified in T0, T1, T2, and T3 group, respectively, while 305, 279, 165 and 298 novel miRNAs were predicted respectively (Table S1). The expression levels of these miRNAs are different, but several miRNAs, such as dre-miR-100-5p, dre-let-7e, dre-miR-101a, dre-miR-143 and dre-miR-146a, were all highly expressed in different levels of hypoxia stress (Table 3). According to the exon, intron and repeat sequence information in genome annotation, and the sRNA information in mirBase and Rfam databases, the small RNA was classified and annotated (Table 4).

Validation of Selected miRNAs by Real-Time PCR
In order to validate the results of high-throughput sequencing, 16 miRNAs were randomly selected from different groups for real-time PCR. The results showed that the expression trends of these miRNAs were basically consistent with the high-throughput sequencing results, indicating that the sequencing results were highly reliable and could be used for subsequent analysis (Figure 2).

Target Gene Prediction and Functional Annotation of Differentially Expressed miRNAs
MiRanda and qTar software were used to predict the target genes of differentially expressed miRNAs in silver carp under hypoxia stress, and the corresponding number of miRNA target genes was 25,146. The number of target genes for known miRNAs ranged from 29 (dre-miR-736) to 996 (dre-miR-17a-3p), while the number of target genes for novel miRNAs ranged from 39 (novel-878) to 1140 (novel-602). Among these target genes, a series of genes related to hypoxia regulation, such as hypoxia-inducible factor-1alpha, vascular endothelial growth factor and hypoxia up-regulated protein, are included.
Next, GO and KEGG enrichment analyses were performed. GO enrichment analysis showed the similar results in the comparison of T0 vs. T1, T0 vs. T2 and T0 vs. T3. The target genes of miRNAs were identified for 55 enriched GO categories in terms of cellular component (CC), molecular functional (MF) and biological process (BP). Of these GO terms, integral component of membrane and membrane part in CC were significantly enriched (p ≤ 0.05). The biological process category of GO terms showed target genes were classified into 26 groups; the top of the most enriched GO term list were cellular process, biological regulation and single-organism process. The cellular component category of GO terms showed target genes were classified into 17 groups, the top of the most enriched GO term list were membrane, macromolecular complex and organelle. In the category of molecular function, binding, catalytic activity and signal transducer activity were the most significantly enriched with respect to 12 groups (Figure 3).

Target Gene Prediction and Functional Annotation of Differentially Expressed miRNAs
MiRanda and qTar software were used to predict the target genes of differentially expressed miRNAs in silver carp under hypoxia stress, and the corresponding number of miRNA target genes was 25,146. The number of target genes for known miRNAs ranged from 29 (dre-miR-736) to 996 (dre-miR-17a-3p), while the number of target genes for novel miRNAs ranged from 39 (novel-878) to 1140 (novel-602). Among these target genes, a series of genes related to hypoxia regulation, such as hypoxia-inducible factor-1alpha, vascular endothelial growth factor and hypoxia up-regulated protein, are included.
Next, GO and KEGG enrichment analyses were performed. GO enrichment analysis showed the similar results in the comparison of T0 vs. T1, T0 vs. T2 and T0 vs. T3. The target genes of miRNAs were identified for 55 enriched GO categories in terms of cellular component (CC), molecular functional (MF) and biological process (BP). Of these GO terms, integral component of membrane and membrane part in CC were significantly enriched (p ≤ 0.05). The biological process category of GO terms showed target genes were classified into 26 groups; the top of the most enriched GO term list were cellular process, biological regulation and single-organism process. The cellular component category of GO terms showed target genes were classified into 17 groups, the top of the most enriched GO term list were membrane, macromolecular complex and organelle. In the category of molecular function, binding, catalytic activity and signal transducer activity were the most significantly enriched with respect to 12 groups (Figure 3). Animals 2021, 11, 2917 9 of 14 According to the results of KEGG enrichment analysis, T0 vs. T1, T0 vs. T2 and T0 vs. T3 also showed the same enrichment trend. The 34 most significantly enriched pathways According to the results of KEGG enrichment analysis, T0 vs. T1, T0 vs. T2 and T0 vs. T3 also showed the same enrichment trend. The 34 most significantly enriched pathways involving in 12,389, 16,446 and 14,052 target genes were divided into five clus-ters, respectively. The top five significantly enriched KEGG pathways include transport and catabolism, cell growth and death, immune system, signal transduction and lipid metabolism (Figure 4, Tables S2-S4). involving in 12389, 16446 and 14052 target genes were divided into five clusters, respectively. The top five significantly enriched KEGG pathways include transport and catabolism, cell growth and death, immune system, signal transduction and lipid metabolism ( Figure 4, Table S2, S3 and S4).

Discussion
Hypoxia is a key factor limiting the survival rate of silver carp with high cultivation density. So far, there have been many studies on the response mechanism of silver carp under hypoxia stress [7]. Since miRNA was first discovered in Caenorhabditis elegans, increasing studies have confirmed that miRNA played an important role in the response of many aquatic organisms to stress [23]. However, the involvement of miRNAs in hypoxia tolerance in silver carp has not been extensively studied [2]. Therefore, we used the Illumina Noveseq high-throughput sequencing platform to study the small RNA in the mixed samples of heart, brain, liver and gill of silver carp under different levels of hypoxia stress, in order to lay a foundation for the further study of the regulatory role of miRNAs in the hypoxic stress response of silver carp.
A total of 229 known miRNAs and 391 putative novel miRNAs were identified, which provided important information for further study of the mechanism of miRNAs in silver carp. By comparing the four hypoxic treatment groups, differentially expressed miRNAs upon hypoxia exposure were identified. Many of these miRNAs have been shown to be closely related to the regulation of hypoxic stress in aquatic organisms. Under hypoxia stress, 14 differentially expressed miRNAs were found in medaka (Oryzias latipes), and miR-204-5p was identified to be able to regulate apoptosis-related genes [24,25]. In this study, miR-204-5p showed an obviously high expression level, which suggested that miR-204-5p played an important role in regulating apoptosis of silver carp under hypoxia stress. The brain, liver and gonad of marine medaka were used as research materials, and highthroughput miRNA group sequencing was carried out [26]. It was found that the expression of let-7a, miR-122 and miR-9-3p was significantly down-regulated in female marine medaka under hypoxia, while miR-2184 was significantly up-regulated [27]. However, so far, there were few studies on the expression of miRNAs in silver carp under hypoxia stress.
In the experiment of cerebral ischemia and hypoxia in rats, let-7f was found to be the upstream regulatory factor of NDRG3. When the expression of let-7f increased, the expression of NDRG3 decreased. Under hypoxic conditions, stable NDRG3 protein could promote angiogenesis and cell growth by activating the Raf-ERK pathway [28]. The experimental results confirmed the new mechanism of let-7f regulating the expression of NDRG3. The miR-125a has also been proven to play an important role in regulating animal-life activities under hypoxia stress [27]. In the study of the effect of hypoxia and reoxygenation on rat cardiomyocytes, miR-125a-5p can regulate the scavenger receptor B1 (Scarb1) gene, and the expression of miR-125a-5p in cardiomyocytes was affected by hypoxia and reoxygenation. When the expression of miR-125a-5p was inhibited, the expression of Scarb1 was significantly up-regulated, which played a protective role in cardiomyocytes [21]. MiR-146a has been confirmed to be involved in the occurrence and progression of a variety of cancers and played a key role in cell apoptosis [19]. The studies on miR-146a under hypoxia stress mainly focus on some human diseases. In the pathogenesis of osteoarthritis, it is confirmed that miR-146a played a key role in the autophagy of chondrocytes under hypoxia. The experimental results showed that miR-146a can promote autophagy by inhibiting the expression of Bcl-2 under hypoxia [29].
Except for the above differentially expressed miRNAs, other miRNAs with high expression levels such as dre-miR-143 and dre-miR-21 also played an important role in regulating the life activities of fish under hypoxia stress [30]. In the physiological adaptation of Atlantic cod (Gadus morhua) and Pacific whiteleg shrimp (Litopenaeus vannamei) under hypoxia stress, it was found that miR-143 could directly target the key enzyme-hexokinase in the mediation of glycolysis pathway, while hexokinase could accelerate glycolysis and produce ATP under hypoxia conditions [31]. In the study of the effect of hypoxiareoxygenation on hepatocyte injury, miR-21 has been found to inhibit the activity of lactate dehydrogenase (LDH), thereby reducing the apoptosis of liver cells, and ultimately reducing the damage of liver cells [32][33][34]. In this study, the high expression of miR-21 was speculated to play an important regulatory role in hypoxia stress response of silver carp.
GO and KEGG enrichment analysis of differentially expressed miRNAs target genes showed that differentially expressed miRNAs target genes were mainly involved in the regulation of hypoxia stress in silver carp through biological regulation, growth and death, immune system and signal transduction. In addition, relevant study has also found that miRNA played an important role in fat metabolism in extreme environments, which wasconsistent with the significant enrichment of lipid metabolism pathways in this study [35]. The results showed that under hypoxia conditions, most of the differentially expressed target genes of miRNAs were related to cell apoptosis, immunity and lipid metabolism [36,37]. This study explores the expression of miRNA in silver carp under different levels of hypoxia stress, and provides the reference for carrying out miRNA-target gene interaction in the later stage, improving the survival ability of silver carp under hypoxia environment.

Conclusions
Overall, the miRNAs expression of silver carp under different hypoxia stress levels was studied by high-throughput sequencing technology, which proved that miRNAs could be involved in the hypoxia stress of silver carp. Many miRNAs involved in the regulation of hypoxia in aquatic organisms, such as miR-204-5p, miR-143 and miR-125a-5p, have been identified. In addition, the target of miRNAs were predicted and functionally by GO and KEGG analysis, suggesting that these miRNAs may play a role in cell apoptosis and immunity. Our results provided new insights to further study the miRNAs regulatory mechanisms and molecular characteristics of hypoxia response in silver carp.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
All authors declare no conflict of interest.