Differential Expression of miRNAs in the Respiratory Tree of the Sea Cucumber Apostichopus japonicus Under Hypoxia Stress

The sea cucumber, an important economic species, has encountered high mortality since 2013 in northern China because of seasonal environmental stress such as hypoxia, high temperature, and low salinity. MicroRNAs (miRNAs) are important in regulating gene expression in marine organisms in response to environmental change. In this study, high-throughput sequencing was used to investigate alterations in miRNA expression in the sea cucumber under different levels of dissolved oxygen (DO). Nine small RNA libraries were constructed from the sea cucumber respiratory trees. A total of 26 differentially expressed miRNAs, including 12 upregulated and 14 downregulated miRNAs, were observed in severe hypoxia (DO 2 mg/L) compared with mild hypoxia (DO 4 mg/L) and normoxic conditions (DO 8 mg/L). Twelve differentially expressed miRNAs were clustered in severe hypoxia. In addition, real-time PCR revealed that 14 randomly selected differentially expressed miRNAs showed significantly increased expressions in severe hypoxia and the expressions of nine miRNAs, including key miRNAs such as Aja-miR-1, Aja-miR-2008, and Aja-miR-184, were consistent with the sequencing results. Moreover, gene ontology and pathway analyses of putative target genes suggest that these miRNAs are important in redox, transport, transcription, and hydrolysis under hypoxia stress. Notably, novel-miR-1, novel-miR-2, and novel-miR-3 were specifically clustered and upregulated in severe hypoxia, which may provide new insights into novel “hypoxamiR” identification. These results will provide a basis for future studies of miRNA regulation and molecular adaptive mechanisms in sea cucumbers under hypoxia stress.

MicroRNAs (miRNAs) are endogenous noncoding small RNAs with 22 nucleotides, which play important roles in various physiological processes, including proliferation, differentiation, apoptosis, and immune responses (Bartel 2004;Chen et al. 2004;Croce and Calin 2005;Pedersen et al. 2007). By binding to target mRNA transcripts, miRNAs can reversibly inhibit translation of mRNAs and/or target them for degradation, thus they are now recognized as key regulators of gene expression (Bartel 2004). MiRNAs are sensitive to environmental changes, and their differential expression is important for the adaption of organisms to the environment (Kim et al. 2014). A previous study showed that miRNAs may play a key role in achieving a hypometabolic state among stress-tolerant animals (Chen et al. 2013). In recent years, a multitude of reports have demonstrated that specific miRNAs are involved in the hypoxic response and contribute to the regulation of hypoxia-related genes, such as HIF, VEGF, and argo-naute2 (Hua et al. 2006;Donker et al. 2007;Camps et al. 2008;Fasanaro et al. 2008;Guimbellot et al. 2009;Yeh et al. 2013). The regulation of organisms to cope with hypoxia stress is under the control of specific hypoxia-inducible miRs, which are also known as "hypoxamiRs" (Chan et al. 2012;Greco et al. 2014). MiR-210, the most accepted master hypoxamir (Chan et al. 2012), and other multiple hypoxamirs, including miR-200b (Ye et al. 2011), miR-199a (Rane et al. 2009), miR-21 (Cheng et al. 2010), miR-23a (Sayed et al. 2010), and miR-181c (Zhang et al. 2012), have been demonstrated to directly target important gene transcripts that regulate cell proliferation, DNA repair, antiinflammatory factors, apoptosis, and angiogenesis, among the adaptions in organisms toward hypoxic conditions. A few miRNA expression profiles have been successfully constructed in aquatic species under hypoxia stress, including river prawn (Macrobrachium nipponense), medaka (Oryzias melastigma), darkbarbel catfish (Pelteobagrus vachelli), and intertidal snail (Littorina littorea) (Biggar et al. 2012;Lau et al. 2014;Sun et al. 2016;Zhang et al. 2016). In A. japonicus, miRNA expression profiles constructed by high-throughput sequencing methods have been applied in different studies, such as regeneration, aestivation, skin ulceration syndrome (SUS), and specific tissues Chen et al. 2013;Wang et al. 2015;Sun et al. 2017). However, "hypoxamiRs" have not yet been reported in sea cucumber and very few have been reported in studies of Echinoderms. More efforts are needed to uncover the regulation mechanisms of miRNAs in response to hypoxia stress in Echinoderms.
In this work, we present for the first time an analysis of the global profile of small RNAs in sea cucumbers using Illumina sequencing technology, and compare them in hypoxia (DO 2 mg/L and DO 4 mg/L) and normoxic (DO 8 mg/L) states. We focus on the respiratory tree because it is the major tissue responsible for respiration and metabolism of sea cucumbers under adverse conditions. The main objectives are to identify and characterize miRNAs, which may play a major role in regulating genes related to hypoxia and stress resistance. Real-time PCR was used to confirm differentially expressed miRNAs. Our findings will be helpful in the further study of sea cucumber biomarkers and hypoxamiR identification under hypoxia stress and will provide important new insights into the molecular mechanisms in sea cucumbers for coping with hypoxic environments.

Animals
Sea cucumbers (body weight 100 6 20 g) were collected from the coast of Weihai, China (East Ocean Science and Technology Co., Ltd). After weighing, they were acclimated in tanks containing aerated sand-filtered seawater (salinity 30‰, pH 8.0) at 15 6 1°for 1 wk before use and were fed once a day at 11:00 AM. Remaining feed was removed daily during the acclimation and experimental periods. All animals were then divided randomly into three groups. One group of sea cucumbers was maintained as the control group (DO8 group) in water with a high DO level of 8 mg/L, and the other two groups were kept in sea water with nitrogen aeration to gradually decrease DO to 4 mg/L (DO4 group) and 2 mg/L (DO2 group) for 3 d. Following exposure, specified quantities of sea cucumbers were dissected promptly and tissues were sampled to be preserved in liquid nitrogen and stored at 280°until subsequent analysis. In this study, three biological replicates were set in each group. The respiratory trees of three sea cucumbers were sampled for RNA extraction in each biological replicate. For instance, DO2_1, DO2_2, and DO2_3 were set in the DO2 group, which means that the nine healthy sea cucumbers were selected for further experimentation after exposure to DO levels of 2 mg/L for 3 d. No dead or dying animals were used in this section, ensuring the changes in miRNA expression were due to hypoxia.
Small RNA library construction and sequencing Respiratory tree tissue samples from nine sea cucumbers (three biological replicates · three sea cucumbers per biological replicate) randomly selected from each group were used in this study. Total RNA from each sample was extracted using an Animal Tissue RNA Purification Kit (LC Sciences, Houston, TX) according to the manufacturer's instructions. Total RNA quality was checked with a Bioanalyzer 2100 (Agilent, Santa Clara, CA) with RNA integrity number .7.0. Total RNA of the sea cucumbers was mixed in equal amounts into nine pooled samples [severe hypoxia stress (DO 2 mg/L): DO2_1, DO2_2, DO2_3; mild hypoxia stress (DO 4 mg/L): DO4_1, DO4_2, DO4_3; normoxic condition (DO 8 mg/L): DO8_1, DO8_2, DO8_3]. Subsequently, a small RNA library was prepared according to the protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA). Next, a 16-30 nt size range of RNA was purified from 15% polyacrylamide gels and then ligated sequentially to 59 and 39 adapters. Reverse transcription was performed followed by polymerase chain reaction (PCR) amplification. The purified PCR products were sequenced by Illumina Hiseq2500 Analysis (LC-BIO, Hangzhou, China).

Sequence data analysis
The raw reads obtained from Hiseq2500 sequencing were trimmed. Clean data were processed by removing low-quality reads, reads with 59 primer contaminants, reads without 39 adapters, reads without insert fragments, reads of ,18 nt, and reads containing poly(A). ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA) and small nucleolar RNA (snoRNA) were identified and removed by blasting against the GenBank database (http://blast. ncbi.nlm.nih.gov) and the Rfram database (http://rfam.xfam.org/). Because no miRNA information on the sea cucumber was in the miRBase21.0, the remaining clean reads were aligned to all known precursor/mature miRNAs of all animal species in miRBase 21.0 with #2 matches. The miRNAs with the highest expression for each mature miRNA family were selected as the temporary miRNA database. Clean data were aligned to the above temporary miRNA database, and the expression of miRNA was determined after summing the count of reads aligned to the temporary miRNA database with #2 mismatches. Finally, the precursor of the identified miRNAs was predicted, and molecules that could not fold into a hairpin structure were regarded as pseudo-miRNA. Potentially novel miRNAs were identified using MIREAP (http://sourceforge.net/projects/mireap/) with stem-loop structure prediction, as previously described (Chen et al. 2013).
Differentially expressed miRNAs were calculated by comparing the miRNA expression between control (DO8 group) and treatment (DO2 and DO4 group) samples. The expression of miRNA in control and treatment samples was normalized to determine the expression of reads per million (RPM). The normalization formula was: Normalized expression = Actual miRNA count/Total count of clean reads · 1,000,000. The final RPM between the biological replications was obtained through averaging. The log 2 (fold-change) and P-value were calculated from the normalized expression. Fold-change = Treatment-mean/ Control-mean and P-value was: If the normalized expression was zero, it was changed to 0.01. A P , 0.05 was used as the criterion to determine the significance of the difference in miRNA expression. A heatmap of the differentially expressed miRNAs was constructed using the pheatmap package in R (version 3.1.3). All unique and shared differentially expressed miRNAs in the three groups were presented in a Venn diagram using Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/index.html).

miRNA target prediction and gene ontology (GO) enrichment analysis
The 39UTRs from the sea cucumber transcriptome assembly (Sun et al. 2011;Du et al. 2012) were extracted as the candidate database to predict the target genes using RNAhybrid software. Two computational target prediction algorithms (TargetScan 50 and miRanda 3.3a) were used to predict the genes targeted by differentially expressed miRNAs. The target genes were removed according to the following criteria: (1) context score percentile ,50 calculated by TargetScan algorithms; and (2) max energy .210 calculated by miRanda algorithms. The intersection of these two softwares was taken as the ultimate genes targeted by differentially expressed miRNAs. The network of predicted genes of miRNAs was shown using Cytoscape 2.8.3 software (Shannon et al. 2003).
The target genes of differentially expressed miRNAs were mapped to GO terms in the database (http://www.geneontology.org) using the program Blast2GO and the Kyoto Encyclopedia of Genes and Genomes Pathway database (KEGG) (http://www.genome.jp/kegg) for GO and KEGG analyses. The GO and pathway terms conforming to a P # 0.05 through Bonferroni's correction were defined as significantly enriched GO terms and pathways.
Quantitative miRNA real-time PCR assay For miRNA analyses, 1 mg of total RNA was reverse-transcribed with miRNA-specific stem-loop RT primers and reverse transcriptase M-MLV (RNase H 2 ) (Takara, Shiga, Japan) according to the experimental protocol. Stem-loop reverse transcription primers were designed following a method described by Chen et al. (2005). The reaction proceeded for 10 min at 25°, 50 min at 42°, followed by 5 min at 85°, and a final hold at 4°. The cDNA was amplified by real-time PCR using Platinum SYBR Green qPCR SuperMix-UDG (Invitrogen, Carlsbad, CA) with miRNA-specific forward and reverse primers (Supplemental Material, Table S1) As an internal control to normalize for technical variations, 5.8 s rRNA was also amplified. All reactions were performed on three biological replicates, each of them being run three times. The 2 2DDCT method was used to analyze the relative expression level of miRNA, and the level of significance was analyzed by SPSS statistics 18.0 software.

Data availability
All raw data are available in the GEO database (GSE100603). Table  S1 contains selected miRNAs and their primer sequences used for real-time PCR. Table S2 includes predicted target annotations of differentially expressed miRNAs. Figure S1 presents length distributions of reads in nine libraries.

Small RNA library construction
To identify differentially expressed miRNAs of sea cucumbers under hypoxia stress, nine small RNA libraries (treatments: DO2_1, DO2_2, and DO2_3; DO4_1, DO4_2, and DO4_3; controls: DO8_1, DO8_2, and DO8_3) were constructed from the respiratory tree of sea cucumber. All raw data have been submitted to the GEO database (GSE100603). High-throughput Illumina Hiseq2500 sequencing of these small RNA libraries yielded a total of 10,215,631 6 689,938 (DO2), 10,901,071 6 748,868 (DO4), and 11,064,987 6 476,564 (DO8) high-quality clean reads, ranging from 16 to 30 nt after trimming the adapter sequences (  (Table  1). After removing the above small RNAs, the other RNAs were further analyzed for the identification of sea cucumber respiratory tree miRNAs against the latest miRBase release version 21.0. A total of 99, 180, and 299 unique miRNAs were searched in the three DO2 libraries; 142, 121, and 185 unique miRNAs were searched in the three DO4 libraries; and 103, 118, and 117 unique miRNAs were searched in the three DO8 libraries. The curves of the corresponding length distributions of reads show similar trends between these libraries by following a typical distribution pattern, and the majority of sequences were 22 nt in length ( Figure S1). A Venn diagram shows the number of specific and common miRNAs in the three groups ( Figure 1). Among them, 75 were classified as common miRNAs, seven as specific miRNAs in the DO2 library, seven as specific miRNAs in the DO4 library, and six as specific miRNAs in the DO8 library.

Different expression profiles of miRNAs and real-time PCR validation
In this study, 187 known miRNAs and 289 novel miRNAs were identified in the sea cucumber respiratory tree. The differentially expressed miRNAs were selected using a t-test with P , 0.05. Valid miRNAs were identified as those detected in all three groups from the same treatment condition, thus a total of 26 nonrepetitive miRNAs were identified to be differentially expressed in the pairwise comparison among the three treatments with significant levels (P , 0.05) ( Table 2). Among the differentially expressed miRNAs, 12 upregulated miRNAs and 14 downregulated miRNAs were observed in severe hypoxia (DO2 group) compared with mild hypoxia (DO4 group) and normoxic conditions (DO8 group). Seven upregulated miRNAs and 12 downregulated miRNAs were obtained in severe hypoxia n compared with normoxic conditions; three upregulated miRNAs and two downregulated miRNAs were obtained in mild hypoxia compared with normoxic conditions; six upregulated miRNAs and five downregulated miRNAs were obtained in the severe hypoxia compared with mild hypoxic conditions (Table 2). In addition, three novel miRNAs (novel-miR-1, novel-miR-2, and novel-miR-3) were identified. The log 2 Ratio log 2

DO4ðmeanÞ
DO8ðmeanÞ of novel-miR-3 were .8, indicating that the three novel miRNAs were significantly upregulated under low DO conditions (2 mg/L and 4 mg/L). Moreover, Aja-miR-7977 and Aja-miR-2835 also meet the criterion of "absolute value of log 2 Ratio .8," indicating that Aja-miR-7977 has significantly decreased expression and Aja-miR-2835 has significantly increased expression in severe hypoxia (DO 2 mg/L) ( Table 2). Thus, they might be important to sea cucumber during hypoxia stress according to the analysis of high-throughput sequencing data. The heatmap of 26 differentially expressed miRNAs (P , 0.05) in three experiment groups are illustrated in Figure 2. As Figure 2 shows, 12 of the miRNAs are clustered in severe hypoxic conditions (DO2 group), including novel-miR-1, novel-miR-2, novel-miR-3, Aja-miR-1, and Aja-miR-2008; six of them are clustered in mild hypoxic conditions (DO4 group), including Aja-miR-71b, Aja-miR-71-5p, and Aja-miR-125-5p; eight of them are clustered in normoxic conditions (DO8 group), including Aja-miR-2011-3p, Aja-miR-31-5p, and Aja-miR-200-3p. We speculated that some of them may be crucial for regulating hypoxia-associated gene expression and play important roles in stress resistance. Notably, the three novel miRNAs (novel-miR-1, novel-miR-2, and novel-miR-3) identified by this study were all clustered in the severe hypoxic condition ( Figure  2), which may provide new insights into novel hypoxia biomarker identification.
GO and pathway enrichment analysis for target genes of the miRNAs To further understand the biological function of miRNAs, the putative target genes of differentially expressed miRNAs were predicted and the annotations are shown in Table S2. Subsequently, GO analysis was used to identify enriched functional groups (P , 0.05) (Figure 4). These target genes predominantly participate in biological processes, cellular component and molecular function. Within the terms of biological processes, the three most frequent categories were "regulation of transcription," "oxidation-reduction processes," and "transport." The three most highly represented molecular component categories were "cytoplasm," "nucleus," and "integral component of membrane." Finally, the three most abundant molecular function categories were "metal ion binding," "ATP binding," and "protein binding." Enriched metabolic pathways and signal transduction pathways were identified and are listed in Figure 5. Twenty-six significantly enriched pathways for target genes (P , 0.05) mainly involved in "Spliceosome," "Regulation of actin cytoskeleton," "Acute myeloid leukemia," "Transcriptional misregulation in cancers," "VEGF signaling pathway," "Lysosome," "GnRH signaling pathway," "Systemic lupus erythematosus," "Endocytosis," and "Tuberculosis" were screened ( Figure 5).

DISCUSSION
The miRNA expression profiles of the respiratory tree in sea cucumber under different DO levels and obtained using high-throughput sequencing technology are reported in this study. A total of nine libraries were constructed, and 90 million reads were obtained. These results will not only play a complementary role in improving the amount of information for the small RNA genome of the sea cucumber but also provide the basis for miRNA regulation in stress resistance and

28.6439
Values indicate the means 6 SD (n = 3). Means not sharing a given letter (a and b) differ significantly (P , 0.05).
hypoxamiR identification. In our study, 19 miRNAs were significantly differentially expressed between the comparison of severe hypoxia and normoxic conditions (DO2 vs. DO8 group); five between mild hypoxia and normoxic conditions (DO4 vs. DO8 group); and 11 between severe hypoxia and mild hypoxia (DO2 vs. DO4 group). The hierarchical clustering result of differentially expressed miRNAs showed that 12 were clustered in the DO2 group, six were clustered in the DO4 group, and eight were clustered in the DO8 group. Real-time PCR was used to verify the expression profile. Fourteen miRNAs were significantly differentially expressed, most of which have been demonstrated to play an important role in oxidative response, apoptosis, cell proliferation, and migration (Nohata et al. 2011;Boulias and Horvitz 2012;Jiang et al. 2016). Combined with the published conclusions in previous studies (Liu et al, 2010a;Liu et al, 2010b;Liu et al, 2015;Nohata et al, 2011;Li et al, 2012), sequencing results, and validation results of expression levels, we speculated that key miRNAs, including Aja-miR-1, Aja-miR-2008, Aja-miR-184, and Aja-miR-31-5p, may play important roles in regulating hypoxia-associated gene expression, thereby enhancing the ability of stress resistance. This result was expected, as a complicated process is involved in the adaption to hypoxia stress.
Three novel miRNAs (novel-miR-1, novel-miR-2, and novel-miR-3) were identified in this study and were specifically clustered in the DO2 group as analyzed by hierarchical clustering (Figure 2). The real-time PCR results showed that novel-miR-1 and novel-miR-3 were both upregulated about twofold and novel-miR-2 was upregulated about threefold under severe hypoxia (DO 2 mg/L) compared with normoxic conditions (DO 8 mg/L) and were consistent with sequencing. Thus, we speculated that they may be important in sea cucumber stress response and adaption under hypoxia and they might be identified as new hypo-xamiRs. The key predicted genes of these three novel miRNAs are  shown in Figure 6. Cathepsin-L (CTSL) was the predicted gene of novel-miR-2 and novel-miR-3. It has been reported that the production of CTSL was induced in fibroblasts and KHT-LP1 cells exposed to hypoxia (Anderson et al. 1989;Cuvier et al. 1997). Furthermore, CTSL was able to broadly affect the immune system and might be involved in increased invasion capacity (Cuvier et al. 1997;Lombadi et al. 2005). In addition, fructose-bisphosphate aldolase A (ALDOA) was also the copredicted gene of novel-miR-2 and novel-miR-3. ALDOA, a type of glycolytic enzyme which could convert glucose to lactate, was identified to be upregulated in hypoxic conditions and could be mediated by hypoxia-inducible factor 1 (HIF-1) (Sørensen et al. 2009;Semenza 2012;Leisz et al. 2015). HIF-1 is reported to play important roles in molecular adaption to hypoxia and is primarily responsible for the subsequent activation of several other hypoxia-responsive genes involved in glycolysis, erythropoiesis, catecholamine metabolism, angiogenesis, transposons and iron metabolism, growth suppressor genes, and several other genes directly linked to the absence of oxygen (Bracken et al. 2003;Egg et al. 2013;Semenza 2014). Thus, we analyzed the potential targeted relationship between novel-miR-2 and HIF-1 in sea cucumber by a manual blast search of the sequence. Results showed that novel-miR-2 was predicted to match 6 bp with the seed sequence of HIF-1a in A. japonicus. Furthermore, ALDOA was associated with hypoxia and the antioxidant stress response and has a positive role in ROS production (Ji et al. 2016). Moreover, Armadillo repeat-containing protein 8 (Armc8), the predicted gene of novel-miR-1, was reported to play an important role in regulating cell migration, proliferation, tissue maintenance, and signal transduction (Jiang et al. 2016). Armc8 is also a key component of the CTLH (C-terminal to lissencephaly type-1-like homology motif) complex of mammalian cells, which has been characterized as an FBPase (fructose-1, 6-bisphosphatase)-degrading complex . In conclusion, the three novel miRNAs may play important roles in the response and adaption of sea cucumbers under hypoxia stress by mediating key genes related to hypoxia such as HIF-1 and genes involved in glycolysis, cell migration, and proliferation.
Aja-miR-2008 was clustered in the DO2 group when analyzed by hierarchical clustering (Figure 2), and it was upregulated 7.6-fold in severe hypoxia compared with normoxic conditions as validated by real-time PCR, and the result was consistent with sequencing. MiR-2008 was previously found to be significantly upregulated in diseased sea cucumber by high-throughput sequencing (Zhang et al. 2015). In previous studies, miR-2008 was significantly upregulated after pathogen infection and was speculated to be involved in a SUS outbreak Zhang et al. 2014a). Under hypoxia stress in this study, several sea cucumbers were observed to be distorted with ulcerated skin, the papillae of the sea cucumbers became white during the early stage of hypoxia stress, and then the body walls dissolved. These symptoms are similar to SUS; therefore the upregulated expression of Aja-miR-2008 was not surprising.
According to the miRNA profile of sea cucumbers in this study, the miR-71 family showed a different expression level under hypoxia stress. The expression of Aja-miR-71-5p was upregulated 1.4-fold under severe hypoxia compared with normoxic conditions as validated by real-time PCR and this tendency was in accordance with sequencing, and the expression levels of Aja-miR-71a and Aja-miR-71b were downregulated.  can increase resistance to heat shock and oxidative stress (Boulias and Horvitz 2012). Moreover, miR-71 has previously been linked with the process of aging and the DNA damage response pathway (De Lencastre et al. 2010). In addition, a study associated miR-71 with a response to environmental stresses and nutrient availability through interactions with different targets of the insulin and PI-3K signaling pathways in Caenorhabditis elegans . Therefore, we speculated that the miR-71 family may play important roles in regulating the hypoxia stress response by targeting genes related to stress resistance.
Aja-miR-31-5p was upregulated 3.2-fold under severe hypoxia (DO2 group) compared with normoxic conditions (DO8 group) by real-time PCR. Former studies showed that miR-31 is involved in several signaling pathways, which can mediate cell proliferation, apoptosis, and DNA mismatch repair. For example, miR-31 could increase cell migration, invasion, and proliferation in an ERK1/2 signaling-dependent manner (Meng et al. 2013). Moreover, miR-31-5p has an important role in radiation responses through regulation of hMLH1 expression, which is one of the core DNA mismatch repair genes that also decreases expression under hypoxia stress (Mihaylova et al. 2003;Kim et al. 2014). MiR-31 was also speculated to be involved in a SUS outbreak in a similar way to miR-2008 . Above all, the research suggested that miR-31-5p may play an important role in hypoxia-related genes such as EGFR, VEGF, and HIF. MiR-31-5p could regulate the signaling pathway downstream of EGFR, which has been shown to induce HIF-1a expression in previous studies (Pore et al. 2006;Igarashi et al. 2015). Furthermore, a previous study suggested that miR-31 contributes to the development of head and neck squamous cell carcinoma by impeding factor inhibiting HIF (FIH) to activate HIF under normoxic conditions (Liu et al. 2010b). Thus, Aja-miR-31-5p could play an important role through the EGFR/HIF-1a/VEGF pathways in stress resistance under hypoxic conditions. Aja-miR-184 was upregulated when analyzed by sequencing and specifically clustered into the DO2 group. It was upregulated 3.7-fold under severe hypoxia when validated by real-time PCR. MiR-184 can cause apoptosis when overexpressed and increase cell numbers when inhibited (Foley et al. 2010). Previous studies have shown that miR-184 overexpression inhibits autophagy and exacerbates oxidative damage . Reports by others provide strong evidence to support the idea that miR-184 is an important modulator of stem cell proliferation and growth (Liu et al. 2010a). Furthermore, miR-184 is reported to be involved in many pathways and some of them are related to the oxidative response. For example, miR-184 negatively modulates Wnt signaling both in vivo and in vitro (Takahashi et al. 2015). MiR-184 could inhibit protein expression in human trabecular meshwork cell cytotoxicity, apoptosis, and extracellular matrix via targeting HIF-1a in vivo, and it can also exhibit angiostatic properties through regulating signaling pathways including Akt, TNF-a, and VEGF (Park et al. 2017;Wang et al. 2017). All the results mentioned above showed that Aja-miR-184 may be crucial to the stress response in sea cucumbers by regulating genes related to autophagy, oxidative damage, and cell proliferation under hypoxia stress. Aja-miR-1 was highly expressed by 19.5-fold under severe hypoxia when analyzed by real-time PCR and the tendency was consistent with high-throughput sequencing results. MiR-1 has been demonstrated to be associated with apoptosis-related genes such as heat shock protein (HSP), and indirectly regulates eNOs (He et al. 2011). A former study indicated that miR-1 regulated Hsp60 expression post-transcriptionally and accelerated cardiomyocyte apoptosis through Hsp60 (Shan et al. 2010). Moreover, apoptosis are mediated by miR-1 (Yu et al. 2008). Restoration of miR-1 in cancer cells inhibits cell proliferation, invasion, and migration (Nohata et al. 2011). Previous results indicate that miR-1 directly regulates the levels of VegfA (Vascular endothelial growth factor A) in muscle, which is regulated by the HIF gene in hypoxia (Stahlhut et al. 2012). MiR-1 levels were significantly increased in response to oxidative stress (Tang et al. 2009). The results reported in previous studies are consistent with the present research.
In addition to the results of high-throughput sequencing, real-time PCR validation and hierarchical clustering, Aja-miR-2835, Aja-miR-153, and Aja-miR-153-3p are also found to be key miRNAs, which may be important in regulating hypoxia-associated gene expression. Although their physiological functions remain unclear, their expression patterns and target genes indicate that these miRNAs are likely to play important roles in hypoxia stress tolerance. Further experiments are needed to elucidate the roles of these molecules in this process.
To obtain insight into the potential function of differentially expressed miRNAs under hypoxia stress, we performed a GO and pathway enrichment analysis of their predicted targets. The significantly enriched GO terms were "regulation of transcription," "oxidation-reduction processes," "transport," "proteolysis," "transcription," and "positive regulation of transcription," which are primarily associated with redox, transport, transcription, and hydrolysis. This can be partially explained by the complex biological process involved in sea cucumbers under hypoxia stress, which includes material transport, oxidative stress, and transcription regulation. By searching KEGG, a total of 26 significantly enriched pathways were observed (P , 0.05), with the top 10 significant pathways being "Spliceosome," "Regulation of actin cytoskeleton," "Acute myeloid leukemia," "Transcriptional misregulation in cancers," "VEGF signaling pathway," "Lysosome," "GnRH signaling pathway," "Systemic lupus erythematosus," "Endocytosis," and "Tuberculosis." Since HIF-1a, a typical hypoxia biomarker reported in other species, is known to be related to hypoxia stress and involved in the "VEGF signaling pathway" (Zhang et al. 2014b), and the differentially expressed miRNA, miR-1, was reported to directly regulate the expression of VEGF (Stahlhut et al. 2012), it was not surprising that this pathway was enriched. In addition, hypoxia stress is supposed to be related to immunity, so some pathways related to immune response were also identified, such as "Lysosome," "Spliceosome," and "Endocytosis" (Figure 5).
In conclusion, this study provides a global view of miRNA changes under hypoxia stress. A total of 26 nonrepetitive miRNAs were significantly differentially expressed in three levels of dissolved oxygen, and these molecules were considered to be responsible for the stress resistance of sea cucumbers. The key miRNAs, analyzed by high-throughput sequencing, real-time PCR validation, and hierarchical clustering, such as Aja-miR-1, Aja-miR-2008, Aja-miR-184, and Aja-miR-31-5p, might be significant in regulating hypoxia-associated gene expression. These miRNAs may play important roles in redox, transport, transcription, hydrolysis, and tolerance ability under hypoxic stress, indicating that sea cucumbers cope with hypoxic stress using a complex biological process, which includes material transport, oxidative stress, and transcription regulation. Moreover, the three novel miRNAs (novel-miR-1, novel-miR-2, and novel-miR-3) identified in the present study were all clustered and significantly upregulated in severe hypoxia, which may provide the basis for hypoxamiR research and novel hypoxia biomarker identification. These results provide an insight into the molecular mechanisms that occur in sea cucumber under hypoxia stress. Figure 6 Predicted target genes of three novel miRNAs (novel-miR-1, novel-miR-2, and novel-miR-3).

ACKNOWLEDGMENTS
This work was supported by National Natural Science Foundation of China-Shandong Joint Fund for Marine Science Research Centers (grant no. U16060404), the Agricultural Seed Project of Shandong Province (2016LZGC032), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA11020703), Qingdao applied basic research program (17-1-1-49-JCH), and by the Taishan Scholars Program (Distinguished Taishan Scholars).