Advertisement

Integrated Application of Multiomics Strategies Provides Insights Into the Environmental Hypoxia Response in Pelteobagrus vachelli Muscle

Open AccessPublished:January 10, 2022DOI:https://doi.org/10.1016/j.mcpro.2022.100196

      Highlights

      • First multiomics analysis of mRNA, miRNA, protein, and metabolite in fishes.
      • Liver and muscle were tissue-specific induced by hypoxia.
      • About 70 genes and 16 miRNAs related to hypoxia adaptation were detected.
      • Hypoxia affects muscle function by mediating energy metabolism via HIF pathway.

      Abstract

      Increasing pressures on aquatic ecosystems because of pollutants, nutrient enrichment, and global warming have severely depleted oxygen concentrations. This sudden and significant lack of oxygen has resulted in persistent increases in fish mortality rates. Revealing the molecular mechanism of fish hypoxia adaptation will help researchers to find markers for hypoxia induced by environmental stress. Here, we used a multiomics approach to identify several hypoxia-associated miRNAs, mRNAs, proteins, and metabolites involved in diverse biological pathways in the muscles of Pelteobagrus vachelli. Our findings revealed significant hypoxia-associated changes in muscles over 4 h of hypoxia exposure and discrete tissue-specific patterns. We have previously reported that P. vachelli livers exhibit increased anaerobic glycolysis, heme synthesis, erythropoiesis, and inhibit apoptosis when exposed to hypoxia for 4 h. However, the opposite was observed in muscles. According to our comprehensive analysis, fishes show an acute response to hypoxia, including activation of catabolic pathways to generate more energy, reduction of biosynthesis to decrease energy consumption, and shifting from aerobic to anaerobic metabolic contributions. Also, we found that hypoxia induced muscle dysfunction by impairing mitochondrial function, activating inflammasomes, and apoptosis. The hypoxia-induced mitochondrial dysfunction enhanced oxidative stress, apoptosis, and further triggered interleukin-1β production via inflammasome activation. In turn, interleukin-1β further impaired mitochondrial function or apoptosis by suppressing downstream mitochondrial biosynthesis–related proteins, thus resulting in a vicious cycle of inflammasome activation and mitochondrial dysfunction. Our findings contribute meaningful insights into the molecular mechanisms of hypoxia, and the methods and study design can be utilized across different fish species.

      Graphical Abstract

      Keywords

      Abbreviations:

      ALAS2 (5-aminolevulinate synthase 2), ANGPTL4 (angiopoietin-related protein 4), ATP5A1/B (ATP synthase subunit alpha/beta, mitochondrial), Blvra (biliverdin reductase A), CA (carbonic anhydrase), COX1/2/4 (cytochrome c oxidase 1/2/4), CV (coefficient of variation), Cyb561a3 (cytochrome b ascorbate-dependent protein 3), CYC1 (cytochrome c1, heme protein, mitochondrial), DEG (differentially expressed gene), DEM (differentially expressed metabolite), DEMI (differentially expressed miRNA), DEP (differentially expressed protein), ECE1 (endothelin-converting enzyme 1), FIH-1 (HIF 1-alpha inhibitor), GLUT1 (solute carrier family 2, facilitated glucose transporter member 1), HADHA (hydroxyacyl-CoA dehydrogenase subunit alpha), HB (hemoglobin), HIF-1 (hypoxia-inducible factor 1), HK1 (hexokinase 1), HMOX2 (heme oxygenase 2), HRG-1A (heme transporter hrg1-A), IL-1β (interleukin 1β), KDR1 (kinase insert domain receptor 1), KEGG (Kyoto Encyclopedia of Genes and Genomes), MAPK (mitogen-associated protein kinase), MDH1 (malate dehydrogenase 1), MS/MS (tandem mass spectrometry), Mt-Cyb (Mt-cytochome b), NLR (NOD-like receptor), NOD (nucleotide oligomerization domain), NOS2 (nitric oxide synthase, inducible), OGDH (2-oxoglutarate dehydrogenase), PCA (principal component analysis), PHD2 (egl nine homolog 1), PHD3 (egl nine homolog 3), PLS-DA (partial least-squares discriminant analysis), PPOX (protoporphyrinogen oxidase), QC (quality control), qRT–PCR (quantitative RT–PCR), ROS (reactive oxygen species), SDHa/b (succinate dehydrogenase a/b), TCA (tricarboxylic acid), TEAB (triethylammonium bicarbonate), TFRC (transferrin receptor protein 1), TMT (tandem mass tag), TPM (transcript per million), UGT2A1 (UDP-glucuronosyltransferase 2A10), urod (uroporphyrinogen decarboxylase), UQCRC1/2 (cytochrome b–c1 complex subunit 1/2, mitochondrial), VEGF (vascular endothelial growth factor), VIP (variable importance in projection)
      Aerobic animals use oxygen to generate energy through oxidative phosphorylation within mitochondria. Currently, aquatic ecosystems are increasingly under pressure from pollutants, nutrient enrichment, and global warming, all of which increase the risk of hypoxia in water (
      • Mandic M.
      • Ramon M.L.
      • Gracey A.Y.
      • Richards J.G.
      Divergent transcriptional patterns are related to differences in hypoxia tolerance between the intertidal and the subtidal sculpins.
      ). Aquatic organisms are usually exposed to different concentrations of dissolved oxygen. Severe and sudden hypoxia can be fatal for hypoxia-sensitive fish. Revealing the molecular mechanism of fish hypoxia adaptation will help us find the hypoxic markers caused by environmental stress (
      • Abdelrahman H.
      • ElHady M.
      • Alcivar-Warren A.
      Aquaculture genomics, genetics and breeding in the United States: Current status, challenges, and priorities for future research.
      ).
      To maintain organism function and homeostasis in a hypoxic environment, aquatic organisms, including fishes, respond by various behavioral and physiological adjustments, for example, decreasing swimming rate, increasing respiration rate, increasing oxygen supply, and reducing oxygen consumption (
      • Sappal R.
      • Fast M.
      • Purcell S.
      • MacDonald N.
      • Stevens D.
      • Kibenge F.
      • Siah A.
      • Kamunde C.
      Copper and hypoxia modulate transcriptional and mitochondrial functional-biochemical responses in warm acclimated rainbow trout (Oncorhynchus mykiss).
      ,
      • Zhong X.P.
      • Dan W.
      • Zhang Y.B.
      • Gui J.F.
      Identification and characterization of hypoxia-induced genes in Carassius auratus blastulae embryonic cells using suppression subtractive hybridization.
      ). These processes depend on the division and cooperation among tissues in fishes, such as the liver, muscle, brain, heart, and gill. Relevant research results show that fish have the tissue-specific expression pattern under hypoxic stimulation, which reflects the different metabolic effects of tissues during hypoxia (
      • Zhu C.D.
      • Wang Z.H.
      • Yan B.
      Strategies for hypoxia adaptation in fish species: A review.
      ). For example, the expression pattern of muscle tissue at the transcription level or the protein level is different from that of brain and liver tissue in response to hypoxia, and there are more downregulated genes in muscle tissue (
      • Gracey A.Y.
      • Troll J.V.
      • Somero G.N.
      Hypoxia-induced gene expression profiling in the euryoxic fish Gillichthys mirabilis.
      ,
      • Lardon I.
      • Eyckmans M.
      • Vu T.N.
      • Laukens K.
      • Boeck G.D.
      • Dommisse R.
      1H-NMR study of the metabolome of a moderately hypoxia-tolerant fish, the common carp (Cyprinus carpio).
      ). The different patterns may also be related to the extent of the hypoxic insult that each tissue experiences. For example, the fish heart receives blood directly from the gills and thus will have a better oxygen supply than other tissues (
      • Everett M.V.
      • Antal C.E.
      • Crawford D.L.
      The effect of short-term hypoxic exposure on metabolic gene expression.
      ,
      • Ju Z.
      • Wells M.C.
      • Heater S.J.
      • Walter R.B.
      Multiple tissue gene expression analyses in Japanese medaka (Oryzias latipes) exposed to hypoxia.
      ).
      Recently, various methods in fish were applied for studies of the hypoxia molecular mechanisms, for example, investigations for transcriptional changes using quantitative RT–PCR (qRT–PCR) (
      • Chen N.
      • Chen L.P.
      • Zhang J.
      • Chen C.
      • Wei X.L.
      • Gul Y.
      • Wang W.M.
      • Wang H.L.
      Molecular characterization and expression analysis of three hypoxia-inducible factor alpha subunits, HIF-1α/2α/3α of the hypoxia-sensitive freshwater species, Chinese sucker.
      ,
      • Wang H.
      • Huang C.
      • Chen N.
      • Zhu K.
      • Chen B.
      • Wang W.
      • Wang H.
      Molecular characterization and mRNA expression of HIF-prolyl hydroxylase-2 (phd2) in hypoxia-sensing pathways from Megalobrama amblycephala.
      ), transcriptome (
      • Beck B.H.
      • Fuller S.A.
      • Li C.
      • Green B.W.
      • Zhao H.
      • Rawles S.D.
      • Webster C.D.
      • Peatman E.
      Hepatic transcriptomic and metabolic responses of hybrid striped bass (Morone saxatilis×Morone chrysops) to acute and chronic hypoxic insult.
      ,
      • Gong D.
      • Xu L.
      • Li W.
      • Shang R.
      • Liu S.
      Comparative analysis of liver transcriptomes associated with hypoxia tolerance in the gynogenetic blunt snout bream.
      ,
      • Yang Y.
      • Fu Q.
      • Wang X.
      • Liu Y.
      • Zeng Q.
      • Li Y.
      • Gao S.
      • Bao L.
      • Liu S.
      • Gao D.
      • Dunham R.
      • Liu Z.
      Comparative transcriptome analysis of the swimbladder reveals expression signatures in response to low oxygen stress in channel catfish, Ictalurus punctatus.
      ), and changes in protein using 2D-difference gel electrophoresis (
      • Dowd W.W.
      • Renshaw G.M.
      • Cech Jr., J.J.
      • Kultz D.
      Compensatory proteome adjustments imply tissue-specific structural and metabolic reorganization following episodic hypoxia or anoxia in the epaulette shark (Hemiscyllium ocellatum).
      ,
      • Smith R.W.
      • Cash P.
      • Ellefsen S.
      • Nilsson G.E.
      Proteomic changes in the crucian carp brain during exposure to anoxia.
      ,
      • Wulff T.
      • Jokumsen A.
      • Hojrup P.
      • Jessen F.
      Time-dependent changes in protein expression in rainbow trout muscle following hypoxia.
      ). Biological phenomena induced by environmental stress are changeable and complex, which is also accompanied by complicated regulation of miRNA/gene/protein expression. Multiple omics are a combination of two or more omics methods and can avoid the incomprehensive conclusions of single omics (
      • Sun S.
      • Guo Z.
      • Fu H.
      • Zhu J.
      • Ge X.
      Integrated metabolomic and transcriptomic analysis of brain energy metabolism in the male Oriental river prawn (Macrobrachium nipponense) in response to hypoxia and reoxygenation.
      ). Unfortunately, the study of the responding mechanism in fish undergoing hypoxia using multiomics is negligible, which vastly hinders the comprehensive understanding of this biological phenomenon.
      Pelteobagrus vachelli (Darkbarbel catfish) has become a popular commercial fish species in Asia because of its relatively high yield and affordable price for consumers (
      • Zhang G.
      • Zhao C.
      • Wang Q.
      • Gu Y.
      • Li Z.
      • Tao P.
      • Chen J.
      • Yin S.
      Identification of HIF-1 signaling pathway in Pelteobagrus vachelli using RNA-seq: Effects of acute hypoxia and reoxygenation on oxygen sensors, respiratory metabolism, and hematology indices.
      ). However, because of its high oxygen threshold and high oxygen consumption rate, this species is only found in the rivers, such as the Liaohe, Huaihe, Yangtze, Xiangjiang, Minjiang, and Pearl. These properties indicate that P. vachelli is not only an important aquaculture variety but also a potential model fish for studying the molecular mechanisms of hypoxia. Recently, the miRNA–mRNA regulatory network has been shown to respond to hypoxia in livers of Megalobrama amblycephala (
      • Sun S.
      • Xuan F.
      • Ge X.
      • Zhu J.
      • Zhang W.
      Dynamic mRNA and miRNA expression analysis in response to hypoxia and reoxygenation in the blunt snout bream (Megalobrama amblycephala).
      ) and P. vachelli (
      • Zhang G.
      • Yin S.
      • Mao J.
      • Liang F.
      • Zhao C.
      • Li P.
      • Zhou G.
      • Chen S.
      • Tang Z.
      Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia.
      ). This suggests that the hypoxia-inducible factor 1 (HIF-1) signaling is well known to be involved in response to hypoxia, which is often regulated and controlled by physiological changes. Also, we conducted isobaric tag for relative and absolute quantitation proteomic analysis in P. vachelli livers to uncover hypoxic responsive proteins involved in diverse biological pathways, for example, peroxisome, glycolysis, peroxisome proliferator–activated receptor signaling, and lipid and amino acid metabolism (
      • Zhang G.
      • Zhang J.
      • Wen X.
      • Zhao C.
      • Zhang H.
      • Li X.
      • Yin S.
      Comparative iTRAQ-based quantitative proteomic analysis of Pelteobagrus vachelli liver under acute hypoxia: Implications in metabolic responses.
      ).
      At present, relative to other tissues, studies of hypoxia on muscle tissue are relatively scarce. Therefore, we characterize transcriptomic, miRNAomic, proteomic, and metabolomic changes of P. vachelli muscles under environmental hypoxia simultaneously (Fig. 1). Our findings aim to offer deeper insights into tissue-specific patterns of expression induced by environmental hypoxia.
      Figure thumbnail gr1
      Fig. 1General workflow of bioinformatics analysis. This was conducted to investigate the molecular mechanism of hypoxia adaptation in Pelteobagrus vachelli in terms of four main aspects (changes in post-transcriptional, transcriptional, protein, and metabolite levels).

      Experimental Procedures

      Experimental Design

      Our experiments were performed following the Guidelines for the Care and Use of Laboratory Animals in China. This study was approved by the Nanjing Normal University Animal Ethics Committee (approval no.: SYXK [Jiangsu] 2015-0028). First, 240 P. vachelli (14 ± 1.25 cm in length and 22 ± 1.78 g in weight) were obtained from Nanjing Fisheries Research Institute. They were acclimated for 3 weeks in aquaria with biofiltered water recirculation system (equipped with cooling and heating functions, the volume of 600 l; flow rate of 5 l/min; 26 ± 1 °C) and fed with artificial compound feed containing more than 42.0% protein (Zhenjiang Jiaji Feed Co Ltd). After stopping feeding for 2 days, these fish were used for challenge experiments. All the experimental methods were the same as the description in our previous article (
      • Zhang G.
      • Yin S.
      • Mao J.
      • Liang F.
      • Zhao C.
      • Li P.
      • Zhou G.
      • Chen S.
      • Tang Z.
      Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia.
      ). In brief, the control fish were dissected, and muscle tissue was taken from the same location, which was 2 cm behind the trailing edge of the gill cover. Each control sample (nos.: M0a, M0b, M0c, M0d, M0e, M0f, M0g, M0h, M0i, and M0j; 6.8 mg/l) was made up of five different individual muscles from a total of 50 control fish. Afterward, the dissolved oxygen in the water was reduced from 6.8 to 0.7 mg/l by bubbling in pure nitrogen for 30 to 35 min. After dissolved oxygen was maintained for 4 h at 0.7 mg/l, 50 treated fish (nos.: M4a, M4b, M4c, M4d, M4e, M4f, M4g, M4h, M4i, and M4j; 0.7 mg/l) were rapidly acquired for muscle dissection. Both the control and treatment had 10 biological replicates, respectively, which were used to conduct metabolome analysis. Meanwhile, six of the same harvested samples (nos.: M0a, M0b, M0c, M4a, M4b, and M4c) were used to conduct analyses of the transcriptome, miRNAome, proteome, and qRT–PCR. All samples were stored at −80 °C until use.

      Differentially Expressed Gene Analysis

      The mRNA-Seq method was the same as the description in our previous article (
      • Zhang G.
      • Li J.
      • Zhang J.
      • Liang X.
      • Zhang X.
      • Wang T.
      • Yin S.
      Integrated analysis of transcriptomic, miRNA and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in heterosis.
      ), including RNA extraction, construction of six complementary DNA libraries, Illumina HiSeq 4000 sequencing, de novo assembly, unigene annotation, and functional classification (accession code: SUB7765255) (
      • Patro R.
      • Duggal G.
      • Love M.I.
      • Irizarry R.A.
      • Kingsford C.
      Salmon provides fast and bias-aware quantification of transcript expression.
      ). The expression level for unigenes was calculated by transcripts per million (TPMs) using Salmon (
      • Mortazavi A.
      • Williams B.A.
      • McCue K.
      • Schaeffer L.
      • Wold B.
      Mapping and quantifying mammalian transcriptomes by RNA-Seq.
      ). Differentially expressed genes (DEGs) were selected with p < 0.05 and fold change >2 using R package edgeR (
      • Robinson M.D.
      • McCarthy D.J.
      • Smyth G.K.
      edgeR: A bioconductor package for differential expression analysis of digital gene expression data.
      ). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEG was the same as the description in our previous article, and so as differentially expressed protein (DEP) and differentially expressed metabolite (DEM) (
      • Zhang G.
      • Li J.
      • Zhang J.
      • Liang X.
      • Zhang X.
      • Wang T.
      • Yin S.
      Integrated analysis of transcriptomic, miRNA and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in heterosis.
      ).
      Pathway enrichment analysis using the corrected p < 0.05 as a threshold to find significantly enriched KEGG terms in the input list of DEGs, comparing them to the whole genome background. The calculation formula of p value was as follows:
      p = 1-i=0m-1(Mi)(NMni)Nn


      N represented the number of KEGG annotated genes in P. vachelli; n represented the number of DEGs in N; M represented the number of particular KEGG annotated genes in a genome; and m represented the number of particular KEGG annotated genes expressed differentially in M. After correction for multiple testing, we chose pathways with a p < 0.05 to represent those significantly enriched in DEGs.

      Differentially Expressed miRNA Analysis

      The miRNA-Seq method was the same as the description in our previous article, including construction of six small RNA libraries, Illumina HiSeq 4000 sequencing, miRNA identification, and prediction of target genes of miRNAs (
      • Zhang G.
      • Li J.
      • Zhang J.
      • Liang X.
      • Zhang X.
      • Wang T.
      • Yin S.
      Integrated analysis of transcriptomic, miRNA and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in heterosis.
      ). Differentially expressed miRNA (DEMI) based on standardized deep-sequencing counts was analyzed by Student's t test (p < 0.05) (
      • Li W.
      • Yang Y.
      • Liu Y.
      • Liu S.
      • Li X.
      • Wang Y.
      • Zhang Y.
      • Tang H.
      • Zhou R.
      • Li K.
      Integrated analysis of mRNA and miRNA expression profiles in livers of Yimeng black pigs with extreme phenotypes for backfat thickness.
      ,
      • Zhao C.
      • Zhang G.
      • Yin S.
      • Li Z.
      • Wang Q.
      • Chen S.
      • Zhou G.
      Integrated analysis of mRNA-seq and miRNA-seq reveals the potential roles of sex-biased miRNA-mRNA pairs in gonad tissue of dark sleeper (Odontobutis potamophila).
      ).

      DEP Analysis

      Tandem mass tag (TMT)-based quantitative proteomic analysis was the same as the description in our previous article, including protein extraction and digestion, protein digestion, TMT labeling, low-pH nano-HPLC–tandem mass spectrometry (MS/MS) analysis, database searching, and protein quantitation (
      • Zhang G.
      • Li J.
      • Zhang J.
      • Liang X.
      • Zhang X.
      • Wang T.
      • Yin S.
      Integrated analysis of transcriptomic, miRNA and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in heterosis.
      ).

      Protein Extraction and Digestion

      Muscle samples (M0a, M0b, M0c, M4a, M4b, and M4c) were first frozen to a dry powder using a vacuum freeze drier. The freeze-dried powder was dissolved with 200 μl triethylammonium bicarbonate (TEAB) dissolution buffer (8 M urea/100 mM TEAB) and centrifuged at 16,000g (centrifugal radius, 0.1 m) for 20 min. The obtained supernatant was subsided by adding four volumes of cold acetone containing 10 mM DTT, and the mixture was kept undisturbed for about 2 h. The resuspended powder was incubated at −20 °C for approximately 2 h. After centrifugation at 16,000g for 20 min at 15 °C, the precipitate was collected and mixed with 800 μl cold acetone containing a solution of 10 mM DTT for 1 h at 56 °C to break the protein disulfide bonds. Following a second round of centrifugation at 16,000g for 20 min at 15 °C, the precipitate was collected, dried, and then stored at −80 °C for later use.
      Total protein concentrations were measured using the Bradford method. For each sample, 100 μg of protein was dissolved to 500 μl in a dissolution buffer and then diluted with 500 μl 50 mM NH4HCO3. The samples were reduced by mixing with 10 mM DTT at 56 °C for 30 min, alkylated by mixing with 50 mM iodoacetamide for 30 min in the dark, and diluted fourfold using 10 mM TEAB.

      Protein Digestion and TMT Labeling

      Next, 100 μg of proteins was digested with 1 μg/μl trypsin (Promega), and the resultant peptide mixture was labeled using chemicals from the TMT10plex reagent set (Fisher Scientific). The proteins from M0a, M0b, M0c, M4a, M4b, and M4c were labeled using TMT reagents 126, 127, 128, 129, 130, and 131, respectively. The labeled samples were combined, desalted using the C18 SPE column (Sep-Pak C18; Waters), and dried in vacuo.
      The peptide mixture was dissolved in buffer A (A: 10 mM ammonium formate in water, pH 10.0, adjusted with ammonium hydroxide) and then fractionated by high pH separation using an Aquity UPLC system (Waters Corporation) connected to a reverse-phase column (BEH C18 column, 2.1 × 150 mm, 1.7 μm, 300 Å; Waters Corporation). The high pH separation was performed using a linear gradient from 0% B to 45% B over a course of 35 min (B: 10 mM ammonium formate in 90% acetonitrile, pH 10.0, adjusted with ammonium hydroxide). The column flow rate was maintained at 250 μl/min, and column temperature was maintained at 45 °C. About 12 fractions were collected, and each fraction was dried in a vacuum concentrator for the next step.

      Low-pH Nano-HPLC–MS/MS Analysis

      The specific steps of low-pH nano-HPLC–MS/MS analysis are as follows: The fractions were resuspended with 32 μl solvent C, respectively (C: water with 0.1% formic acid; D: acetonitrile with 0.1% formic acid), separated by nanoLC, and analyzed by online electrospray MS/MS. The experiments were performed on an EASY-nLC 1000 system (Thermo Fisher Scientific) connected to an Orbitrap Fusion mass spectrometer (Thermo Fisher Scientific) equipped with an online nanoelectrospray ion source. About 4 μl of the peptide sample was loaded onto the trap column (Thermo Scientific Acclaim PepMap C18 column; 100 μm × 2 cm), with a flow of 10 μl/min for 3 min and subsequently separated on the analytical column (Acclaim PepMap C18; 75 μm × 25 cm) with a linear gradient, from 5% D to 30% D in 95 min. The column was re-equilibrated at initial conditions for 15 min. The column flow rate was maintained at 300 nl/min, and column temperature was maintained at 45 °C. The electrospray voltage of 2.0 kV versus the inlet of the mass spectrometer was used. The Orbitrap Fusion mass spectrometer was operated in the data-dependent mode to switch automatically between MS and MS/MS acquisition. Survey full-scan MS spectra (m/z 400–1600) were acquired in Orbitrap with a mass resolution of 60,000 at m/z 200. The automatic gain control target was set to 500,000, and the maximum injection time was 50 ms. MS/MS acquisition was performed in Orbitrap with 3 s cycle time, and the resolution was 15,000 at m/z 200. The intensity threshold was 50,000, and the maximum injection time was 150 ms. The automatic gain control target was set to 150,000, and the isolation window was m/z 2. Ions with charge states 2+, 3+, and 4+ were sequentially fragmented by higher energy collisional dissociation with a normalized collision energy of 37%. In all cases, one microscan was recorded using dynamic exclusion of 30 s. MS/MS fixed first mass was set at 110.

      Database Searching

      Database searching was conducted using the present RNA-Seq results (accession code: SUB7765255) by Mascot (version 2.3; Matrix Science). Tandem mass spectra were extracted by Proteome Discoverer software (version 1.4.0.288; Thermo Fisher Scientific). The parameters are set as follows: the mass range of the precursor ion is 350 to 5000 Da; the minimum number of peaks in MS/MS is 10; and the signal-to-noise ratio threshold is 1.5. Charge state deconvolution and deisotoping were not performed. All MS/MS samples were analyzed using Mascot. Mascot was set up to search the UniProt database (taxonomy: Siluroidei, 18,581 entries) assuming the digestion enzyme trypsin (Promega). Maximum missed cleavage is set to 2. Mascot was searched with a fragment ion mass tolerance of 0.050 Da and a parent ion tolerance of 10.0 ppm. Carbamidomethyl of cysteine and TMT6plex of lysine and the N terminus were specified in Mascot as fixed modifications. Oxidation of methionine was specified in Mascot as a variable modification.

      Quantitative Data Analysis

      Use the percolator algorithm to control peptide level false discovery rate lower than 1%. Only unique peptides were used for protein quantification, the method of normalization on protein median was used to correct experimental bias, and the minimum number of proteins that must be observed to allow was set to 1000. A protein was required to contain at least two unique peptides. Nine protein ratios were produced by comparing three M4 and three M0 (i.e., M4a/M0a, M4a/M0b, M4a/M0c, M4b/M0a, M4b/M0b, M4b/M0c, M4c/M0a, M4c/M0b, M4c/M0c). Then the p values were determined using ANOVA and the nine ratios. The mean comparisons of pairs among the nine protein ratios calculated fold changes. Finally, proteins with p < 0.05 and fold change >1.2 were considered to be DEPs (
      • Zhang G.
      • Li J.
      • Zhang J.
      • Liang X.
      • Zhang X.
      • Wang T.
      • Yin S.
      Integrated analysis of transcriptomic, miRNA and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in heterosis.
      ,
      • Chen Z.
      • Wen B.
      • Wang Q.
      • Tong W.
      • Guo J.
      • Bai X.
      • Zhao J.
      • Sun Y.
      • Tang Q.
      • Lin Z.
      • Lin L.
      • Liu S.
      Quantitative proteomics reveals the temperature-dependent proteins encoded by a series of cluster genes in thermoanaerobacter tengcongensis.
      ).

      DEM Analysis

      Metabolome method was the same as the description in a previous study, including 20 metabolite extracts (no.: M0a, M0b, M0c, M0d, M0e, M0f, M0g, M0h, M0i, M0j, M4a, M4b, M4c, M4d, M4e, M4f, M4g, M4h, M4i, and M4j), LC–MS conditions for analysis, and metabolomics data processing method (
      • Cao M.
      • Li C.
      • Liu Y.
      • Cai K.
      • Zhou X.
      Assessing urinary metabolomics in giant pandas using chromatography/mass spectrometry: Pregnancy-related changes in the metabolome.
      ). Data analysis methods include partial least-squares discriminant analysis (PLS-DA) and principal component analysis (PCA). Quantitative identification of metabolite screening and supervised PLS-DA was evaluated by MetaX software (BGI; http://metax.genomics.cn) (
      • Cao M.
      • Li C.
      • Liu Y.
      • Cai K.
      • Zhou X.
      Assessing urinary metabolomics in giant pandas using chromatography/mass spectrometry: Pregnancy-related changes in the metabolome.
      ). The specific steps are as follows: the quantitative information of metabolites comes from the peak area of the first chromatogram of the substance. We use XCMS software (Scripps Research Institute) to extract the intensity information of each substance in each sample and then use MetaX software for quality control (QC) of the extracted data. First, remove low-quality peaks (more than 50% of the QC sample is missing, or the actual sample is missing, and more than 80% of the ions are removed) and then use the K-nearest neighbors method to fill in the missing values and then use probabilistic quotient normalization. Filter the corrected data, that is, filter out ions with a coefficient of variation (CV) >50% in all QC samples (the ions with a CV >50% fluctuate greatly during the experiment, and quantitative analysis of differences is not performed). The variable importance in projection (VIP) was calculated, and metabolites with VIP >1 for both multidimensional statistical analysis and p < 0.05 (Student's t test) were selected as metabolites with significant differences (
      • Wen X.
      • Hu Y.
      • Zhang X.
      • Wei X.
      • Yin S.
      Integrated application of multi-omics provides insights into cold stress responses in pufferfish Takifugu fasciatus.
      ).

      Multiomics Analysis of Transcriptome, miRNAome, Proteome, and Metabolome

      ACGT101-CORR 1.1 (LC Sciences) was used to construct the negative miRNA–mRNA regulatory network (up–down and down–up). In brief, the transcriptome and miRNAome were integrated by combining the expression profiles of DEG and DEMI with miRNA-targeting information. Afterward, the protein information added miRNA–mRNA pairs to the miRNA–mRNA–protein pairs (up–down–down and down–up–up) directly. Second, 145 miRNA–mRNA–protein pairs and 33 DEMs were combined to find mutual KEGG pathways. Finally, Cytoscape software (Institute for Systems Biology) was used to generate an interaction network among the relationships of DEMI, DEG/DEP, and DEM (
      • Wen X.
      • Hu Y.
      • Zhang X.
      • Wei X.
      • Yin S.
      Integrated application of multi-omics provides insights into cold stress responses in pufferfish Takifugu fasciatus.
      ).

      qRT–PCR

      The expression profiles of all 16 DEMIs among the multiomics interaction network and the selected 70 DEGs representing biological processes induced by hypoxia were further analyzed using qRT–PCR. qRT–PCR methods were the same as the description in our previous article (
      • Zhang G.
      • Li J.
      • Zhang J.
      • Liang X.
      • Zhang X.
      • Wang T.
      • Yin S.
      Integrated analysis of transcriptomic, miRNA and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in heterosis.
      ). All primers are shown in supplemental Table S1. 2−△△CT method was used to calculate the expression levels, and the data were subjected to statistical significant analysis (p < 0.05) using SPSS 22.0 (IBM Corp) (Student's t test).

      Results

      Transcriptome Summary

      A total of 28,696 unique genes were identified from six complementary DNA libraries (M0a, M0b, M0c, M4a, M4b, and M4c) representing P. vachelli muscles (supplemental Table S2). The QC analysis of the transcriptome is shown in supplemental Table S3, including sequencing QC overview, trinity assembly results overview, gene/transcript length overview, and gene function checklist. We found that 3511 genes were significantly upregulated, whereas 1033 genes were significantly downregulated in response to hypoxia (p < 0.05 and fold change >2). By performing the KEGG pathway analyses for 4544 DEGs, a total of 22 pathways were identified that showed significant changes after hypoxia (p < 0.05) (Fig. 2A).
      Figure thumbnail gr2
      Fig. 2KEGG pathway enrichment. Enrichment of 4544 differentially expressed genes (DEGs; A), 584 differentially expressed proteins (DEPs; B), and 33 differentially expressed metabolites (DEMs; C) (22, 14, and 14 pathways; p < 0.05). KEGG, Kyoto Encylopedia of Genes and Genomes.
      Among these pathways (Fig. 2A), the “vascular endothelial growth factor (VEGF) signaling pathway” involved in mediating vasopermeability, angiogenesis, and vascular tone was significantly enriched (p < 0.05). However, DEG of “African trypanosomiasis” reflected the biological process of anemia (hemoglobin synthesis, downregulated). Afterward, catabolism and anabolism were the two most represented subclasses. Besides, other enriched subclasses include human immune response and inflammation. Finally, we found that dilated cardiomyopathy and hypertrophic cardiomyopathy represented the impaired muscle function, which was also significantly enriched (p < 0.05). Interestingly, the number of upregulated genes was greater than the number of downregulated genes in the biological processes of vessel generation, catabolism, immune response, inflammation, and muscle dysfunction, whereas the number of downregulated genes exceeded the number of upregulated genes in the processes of anemia and anabolism.

      Proteome Summary

      All MS/MS spectra were processed using Mascot from six protein samples (M0a, M0b, M0c, M4a, M4b, and M4c). TMT analysis of P. vachelli muscle proteome showed 14,943 unique peptides from 41,278 unique spectra in this database (334,625 total spectra) and resulted in 1551 unique proteins. All peptide/protein identifications and quantifications are shown in supplemental Table S2. The QC analysis of the proteome is shown in supplemental Table S3, including mass delta, the distribution diagram of peptide number, peptide length, protein cleavages, and protein mass. All distributions of CV for proteome showed that a total of 97.36% of the identified proteins displayed a ratio percent CV <20%, and the proteome project has good repeatability (supplemental Table S3).
      A total of 584 DEPs (fold change >1.2 and p < 0.05) were reliably quantified using TMT analysis, including nine upregulated and 575 downregulated proteins. KEGG enrichment analysis (p < 0.05) displayed a total of 14 significantly enriched pathways under hypoxia (Fig. 2B). Muscle function, amino acid metabolism, and carbohydrate metabolism were the three most represented subclasses. Also “Alzheimer's disease,” “Huntington's disease,” and “Parkinson's disease” containing mitochondrial-related proteins were significantly enriched by down DEP. Of note, DEP of carbohydrate metabolism referred to aerobic metabolism mainly centered on “pyruvate metabolism” (e.g., citrate synthase, isocitrate dehydrogenase, dihydrolipoyl dehydrogenase, 2-oxoglutarate dehydrogenase [OGDH], succinate dehydrogenase a/b [SDHa/b], malate dehydrogenase 1 [MDH1], downregulated) and the DEP of amino acid metabolism and human disease involved in the amino acid synthesis and muscle dysfunction.

      Correlation of miRNA–mRNA–Protein Regulatory Network

      About 348 conserved miRNAs belonging to 292 miRNA seed-based families were identified in the six small RNA libraries (M0a, M0b, M0c, M4a, M4b, and M4c). The QC analysis of the miRNAome is shown in supplemental Table S3, including an overview of reads from raw data to cleaned sequences, a summary of known and predicted miRNA, the length distribution of sequences, and conservation of the identified miRNA with other species. A total of 39 DEMIs were found (p < 0.05), with 17 upregulated and 22 downregulated miRNAs (supplemental Table S2).
      In our results, according to the DEMI, DEG, and DEP datasets and miRNA-targeting information, 1809 negative miRNA–mRNA regulatory networks with the expression pattern of up–down and down–up were detected (supplemental Table S4). Furthermore, 12 of 20 significant enriched KEGG pathways for 1809 pairs were found to be overlapping for those of transcriptome, also involved in vessel formation, metabolism, immune response, and muscle dysfunction. In addition, similar to the DEG, the numbers of upregulated genes were greater than that of downregulated ones involved in vessel formation, catabolism, immune response, muscle dysfunction, and vice versa for anabolism (Fig. 3A). Finally, we identified 145 negatively correlated miRNA–mRNA–protein triplets (up–down–down and down–up–up) as the key components for our data analysis, including 21 DEMIs and 60 co-DEGs/DEPs in total (Fig. 4A). Moreover, five of 11 significant enriched KEGG pathways for 145 pairs were found to be overlapping for those of the proteome, such as “citrate cycle (tricarboxylic acid [TCA] cycle),” “oxidative phosphorylation,” “tryptophan metabolism,” “glyoxylate and dicarboxylate metabolism,” and “galactose metabolism,” with carbohydrate metabolism as the most commonly represented subclasses. Also, biosynthesis related to “butirosin and neomycin biosynthesis” and “protein export” was significantly enriched by down co-DEG/DEP (Fig. 3B).
      Figure thumbnail gr3
      Fig. 3KEGG pathway enrichment analyses. Analyses of 1809 negative miRNA–mRNA pairs (A) and 145 negative miRNA–mRNA–protein pairs (B) (20 and 12 pathways; p < 0.05). KEGG, Kyoto Encylopedia of Genes and Genomes.
      Figure thumbnail gr4
      Fig. 4About 145 negative miRNA–mRNA–protein pairs. Includes 21 DEMI and 60 co-DEG/DEP (A) and multiomics interaction network including 16 DEMIs, 23 co-DEGs/DEPs, and 17 DEMs (B). DEG, differentially expressed gene; DEMI, differentially expressed miRNA; DEP, differentially expressed protein.

      Metabolome Summary

      A total of 9433 metabolic ion peaks were extracted from P. vachelli muscles in the control group (M0a, M0b, M0c, M0d, M0e, M0f, M0g, M0h, M0i, and M0j) and those in the experimental group (M4a, M4b, M4c, M4d, M4e, M4f, M4g, M4h, M4i, and M4j), including 5152 positive and 4281 negative ion peaks (supplemental Table S2). In addition, we found that 33 unique DEMs (with MS2; 14 upregulated and 19 downregulated) were among the 18 positive and 24 negative DEMs (VIP >1 and p < 0.05; supplemental Table S2).
      The PCA model showed that the positive and negative ion modes were closely clustered. The PLS-DA model of the example set was also established. The model evaluation parameters were positive ion mode (R2 = 0.99 cum, Q2 = 0.78 cum) and negative ion mode (R2 = 0.99 cum, Q2 = 0.87 cum). Aforementioned PCA and PLS-DA results indicated that the model was reliable and stable (supplemental Table S3). Similar to the KEGG analysis of DEG, 14 significant pathways for the biological process of catabolism, anabolism, and muscle dysfunction were enriched by 33 DEMs (p < 0.05). Significant decreases in intermediate metabolites were found in catabolism; however, the upregulated intermediate metabolites were big part of anabolic pathways (Fig. 2C).

      Multiomics Identification of Key miRNAs, Genes, Proteins, and Metabolites

      The hierarchical clusterings of the DEP, DEMI, DEG, and DEM based on the fold-change Z score of M4a/(M4b, M4c, M0a, M0b, M0c), six samples' log10 (norm), six samples' Z score of TPM, and 20 samples' Z score of norm were showed in supplemental Fig. S1. These clusterings indicated that all samples (M4 versus M0) can be divided into two different groups according to control and hypoxia treatment. In general, environmental hypoxia had a significant effect on the overall gene/miRNA/protein/metabolite expression profile of P. vachelli muscle.
      In this study, the transcriptome, miRNAome, proteome, and metabolome of P. vachelli muscle under hypoxia stress were integrated by investigating whether changes in “gene/protein” levels correlated with changes in the corresponding “negative miRNAs” and coupled with mutual KEGG pathways for DEM. This was based on 145 negative miRNA–mRNA–protein pairs, and 33 unique DEMs, 16 DEMIs, 23 co-DEGs/DEPs, and 17 DEMs were involved in the 19 mutual KEGG pathways (Fig. 4B).
      These mutual pathways contained carbohydrate metabolism (i.e., TCA cycle, pentose phosphate pathway, starch, sucrose, pyruvate, butanoate, glyoxylate, and dicarboxylate metabolism), amino acid metabolism (i.e., valine, leucine, isoleucine, cysteine, methionine, alanine, aspartate, glutamate, d-glutamine, and d-glutamate metabolism), other metabolic processes (i.e., purine and nitrogen metabolism, biosynthesis of aminoacyl-tRNA, unsaturated fatty acids, butirosin, and neomycin), membrane transport (i.e., ABC transporters), signal transduction (i.e., mammalian target of rapamycin signaling pathway), and cellular community (i.e., gap junction). In our case, downregulated miR-193b/457a/15b-OGDH, miR-210/338-SDHb, and miR-457a/214-MDH1 led to suppression of aerobic metabolism (e.g., “pyruvate metabolism” and “TCA cycle”), thus making intermediate products rise (e.g., citric acid and S-lactoylglutathione). “Aminoacyl-tRNA biosynthesis” was regulated by downregulated miR-15b/338/457a-Aars/mars and miR-133b/18a/PC-7236-hars to accumulate intermediates (proline, isoleucine, and methionine).

      qRT–PCR Analysis

      To further analyze the hypoxia-induced biological processes in P. vachelli muscles, we performed qRT–PCR analysis on 16 selected DEMIs and 70 DEGs. These representative miRNAs/genes were artificially selected to determine their potential role in environmental hypoxia for P. vachelli muscles. For example, all 16 DEMIs (miR-15b/c, 133b, 181a, 193b, 338, 457a, 1-2, 18a, 30a/d, 210, 214, PC-7236, and let-7b) were from the multiomics identification (Fig. 5), and some DEGs were related to important biological processes (i.e., oxygen sensor, angiogenesis, heme and erythropoiesis, glycolysis, TCA cycle, lipolysis, oxidative stress, mitochondrial dysfunction, nucleotide oligomerization domain [NOD]-like receptor, NF-κB, mitogen-associated protein kinase [MAPK] signaling pathway, cell cycle, and apoptosis). The results of qRT–PCR revealed similar expression patterns between most of these mRNAs/miRNAs and those from transcriptome/miRNAome data (TPM/reads-based expression values) (Table 1 and Figs. 5 and 6). Of these, miR-15c showed minimal expression in control groups (M0a, M0b, and M0c) in the qRT–PCR analysis and was induced by hypoxia stress, which was consistent with the present miRNA-Seq data. Although there were some quantitative differences between the two analysis platforms, our results showed similar trends between the qRT–PCR and RNA-Seq data.
      Figure thumbnail gr5
      Fig. 5Relative miRNA expression of all 16 DEMIs among the multiomics identification for comparison of the M4 versus M0 groups, with respect to quantitative RT (qRT)–PCR and miRNA-Seq. Relative miRNA level (qRT–PCR) is on the left, and relative miRNA level (norm) is on the right of each bar chart. DEMI, differentially expressed miRNA.
      Table 1Relative mRNA expression of 48 selected DEGs for comparison of the M4 versus M0 groups, with respect to mRNA-Seq and qRT–PCR
      AnnotationmRNA-Seq (log2 FC)RegulationqRT–PCR (log2 FC)AnnotationmRNA-Seq (log2 FC)RegulationqRT–PCR (log2 FC)
      Oxygen sensorMitochondrial dysfunction
      FIH2.26Up1.53
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      PPRC1−1.22Down−0.98
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      PHD22.17Up2.13
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      mt-co1−3.87Down−4.41
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      PHD33.45Up3.11
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      mt-co2−3.99Down−3.20
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      AngiogenesisUQCRC2−1.35Down−1.71
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      VEGF2.49Up2.76
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      CYC1−1.24Down−1.98
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      ANGPTL41.12Up2.00
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      mt-Cyb−4.21Down−4.04
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Heme and erythropoiesisAtp5a1−1.27Down−1.14
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      ALAS2−1.99Down−2.46
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      NLR signaling pathway
      TF−4.79Down−3.07
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      NLRC31.36Up2.60
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      TFRC−1.14Down−1.56
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      NLRP31.80Up1.68
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      CA−1.51Down−1.10
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      NLRP121.05Up1.45
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      HBα−1.56Down−2.23
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      IL-1β3.86Up3.40
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      HBβ−1.53Down−2.80
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      NF-κB signaling pathway
      GlycolysisIFNGR11.25Up1.32
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      GLUT11.32Up1.24
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      IKBKG1.53Up2.16
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      PK−1.05Down−1.76
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      NFKBIB1.25Up1.86
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      HK1−1.29Down−1.37
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      CHUK/IKK1.33Up2.14
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      TCA cycleMAPK signaling pathway
      PDK1.88Up1.07
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      MAPKAPK51.43Up2.23
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      PDHA−1.03Down−1.36
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Map4k21.82Up1.24
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      OGDH−1.27Down−0.97
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Map2k41.17Up1.75
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      SDHb−1.16Down−1.97
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Cell cycle
      SDHa−1.28Down−1.56
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      CDK10−1.12Down−0.99
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      MDH1−1.06Down−1.50
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      CCNI−1.23Down−2.19
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      LipolysisCCNG21.89Up1.21
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      PPARα1.21Up1.51
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      DIDO11.94Up1.90
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      PPARβ2.38Up2.97
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      MDC11.87Up1.02
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Oxidative stressCDKN1B2.50Up2.31
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Mn-SOD1.33Up1.01
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      GADD45A2.53Up1.25
      Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Abbreviations: CCNG2, cyclin G2; CCNI, cyclin I; CDK10, cyclin-dependent kinase 10; CDKN1B, cyclin-dependent kinase inhibitor 1B; CHUK, inhibitor of nuclear factor kappa-B kinase subunit alpha; DIDO1, death-inducer obliterator 1; FC, fold change; GLUT1, solute carrier family 2, facilitated glucose transporter member 1; IFNGR1, interferon gamma receptor 1; IKBKG, NF-κB essential modulator; Map2k4, mitogen-activated protein kinase kinase 4; Map4k2, mitogen-activated protein kinase kinase kinase kinase 2; MAPKAPK5, MAP kinase-activated protein kinase 5; Mn-SOD, manganese superoxide dismutase; NFKBIB, NF-κB inhibitor beta; NLRC3, NLR family CARD domain-containing protein 3; PDHA, pyruvate dehydrogenase; PDK, pyruvate dehydrogenase kinase isoform; PK, pyruvate kinase; PPAR, peroxisome proliferator–activated receptor; PPRC1, peroxisome proliferator–activated receptor gamma coactivator-related protein 1; TF, transferrin.
      a Indicates the statistical significance of differential gene expression with p < 0.05 (t test). FC = M4 group (mean)/M0 group (mean).
      Figure thumbnail gr6
      Fig. 6DEG/DEP enrichment of apoptotic pathways including extrinsic, mitochondrial, endoplasmic reticulum stress, and relative mRNA expression of 22 selected DEGs for apoptosis, with respect to mRNA-Seq and quantitative RT–PCR. For abbreviations and explanations, please see the text. DEG, differentially expressed gene; DEP, differentially expressed protein.

      Discussion

      We constructed 145 negative miRNA–mRNA–protein pairs in P. vachelli muscle induced by hypoxia (supplemental Table S4) and later discovered that 16 DEMIs, 23 DEGs/DEPs, and 17 DEMs were involved in 19 mutual KEGG pathways (supplemental Table S5). The predicted miRNA–mRNA–protein–metabolite regulatory networks were more complex than previously thought. A single miRNA could regulate multitarget mRNAs/proteins/metabolites and vice versa (Fig. 4). In performing KEGG pathway analyses for DEG/DEP and multiomics pairs, “signal transduction,” “metabolism,” and “muscle function” were found to be the three most frequently represented subclasses (Figs. 2 and 3). With regard to the fish hypoxia adaptation strategy, we highlighted our particular research using pathway analysis related to the functional clusters, including HIF-1 signaling pathway, energy metabolism, and muscle function. Some DEGs and their regulating miRNAs among the aforementioned clusters representing hypoxic markers were further analyzed using qRT–PCR (Table 1 and Figs. 5 and 6).

      HIF-1 Signaling Pathway

      “HIF-1 signaling” is a major regulator of the hypoxia signaling pathway, leading to similar biochemical and physical reactions, including oxygen sensing, increasing oxygen delivery, and reducing oxygen consumption (
      • Zhu C.D.
      • Wang Z.H.
      • Yan B.
      Strategies for hypoxia adaptation in fish species: A review.
      ).

      Oxygen Sensor

      In recent years, the relationship between several oxygen sensors and hypoxia in some fishes has been elucidated (
      • Wang H.
      • Huang C.
      • Chen N.
      • Zhu K.
      • Chen B.
      • Wang W.
      • Wang H.
      Molecular characterization and mRNA expression of HIF-prolyl hydroxylase-2 (phd2) in hypoxia-sensing pathways from Megalobrama amblycephala.
      ,
      • Geng X.
      • Feng J.
      • Liu S.
      • Wang Y.
      • Arias C.
      • Liu Z.
      Transcriptional regulation of hypoxia inducible factors alpha (HIF-alpha) and their inhibiting factor (FIH-1) of channel catfish (Ictalurus punctatus) under hypoxia.
      ). We found that the transcriptional expressions of the oxygen sensors (HIF 1-alpha inhibitor [FIH-1], egl nine homolog 1 [PHD2], and egl nine homolog 3 [PHD3]) were significantly upregulated, and their regulating miRNAs were significantly downregulated (Table 1), which indicates that oxygen sensors and their regulation of miRNAs play a vital role in response to hypoxia. Similar results were found in P. vachelli liver and brain (
      • Zhang G.
      • Zhao C.
      • Wang Q.
      • Gu Y.
      • Li Z.
      • Tao P.
      • Chen J.
      • Yin S.
      Identification of HIF-1 signaling pathway in Pelteobagrus vachelli using RNA-seq: Effects of acute hypoxia and reoxygenation on oxygen sensors, respiratory metabolism, and hematology indices.
      ).

      Increased Oxygen Delivery

      The pattern of gene upregulation (19 of 23) in the “VEGF signaling pathway” may be involved in promoting their vasopermeability, angiogenesis (VEGF, kinase insert domain receptor 1 [KDR1], angiopoietin-related protein 4 [ANGPTL4], upregulated), and mediating the vascular tone nitric oxide synthase, inducible [NOS2], endothelin-converting enzyme 1 [ECE1], upregulated) (Table 1). Similar activation of the VEGF or angiopoietin-related protein 4 has been shown in P. vachelli liver and Astronotus ocellatus (
      • Zhang G.
      • Yin S.
      • Mao J.
      • Liang F.
      • Zhao C.
      • Li P.
      • Zhou G.
      • Chen S.
      • Tang Z.
      Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia.
      ,
      • Baptista R.B.
      • Souza-Castro N.
      • Almeida-Val V.M.
      Acute hypoxia up-regulates HIF-1alpha and VEGF mRNA levels in Amazon hypoxia-tolerant Oscar (Astronotus ocellatus).
      ).
      It has been widely accepted that oxygen regulation in fish involves improving oxygen uptake and transport to tissues in a hypoxic environment. In P. vachelli, similarly raised blood indices (red blood cell, hemoglobin [HB], and SI/TIBC) and liver mRNA levels (transferrin receptor protein 1 [TFRC], carbonic anhydrase [CA], and erythropoietin) have been reported (
      • Zhang G.
      • Zhao C.
      • Wang Q.
      • Gu Y.
      • Li Z.
      • Tao P.
      • Chen J.
      • Yin S.
      Identification of HIF-1 signaling pathway in Pelteobagrus vachelli using RNA-seq: Effects of acute hypoxia and reoxygenation on oxygen sensors, respiratory metabolism, and hematology indices.
      ). However, in this study, “African trypanosomiasis” was enriched by a downregulated anemia-related gene. We found that the genes or proteins involved in heme synthesis (5-aminolevulinate synthase 2 [ALAS2], heme transporter hrg1-A [HRG-1A], uroporphyrinogen decarboxylase [urod], biliverdin reductase A [Blvra], UDP-glucuronosyltransferase 2A1 [UGT2A1], protoporphyrinogen oxidase [PPOX], and heme oxygenase 2 [HMOX2]) and erythropoiesis (transferrin, TFRC, CA, HBα, and PC-7236-HBβ) were significantly reduced (Table 1 and Fig. 5). This may be because the muscle tissue is not the hematopoietic organ of the fish. However, among the more highly upregulated proteins in the Danio rerio muscle during hypoxia were two variants of HBα, which may be due to different period and level of hypoxia or various regulatory mechanisms in different fish muscle (
      • Chen K.
      • Cole R.B.
      • Rees B.B.
      Hypoxia-induced changes in the zebrafish (Danio rerio) skeletal muscle proteome.
      ).

      Reduce Oxygen Consumption

      In this study, upregulation of solute carrier family 2, facilitated glucose transporter member 1 (HIF-1α target gene) suggested a promotion in anaerobic metabolism in muscle, similar to previously reported results in P. vachelli liver (
      • Zhang G.
      • Yin S.
      • Mao J.
      • Liang F.
      • Zhao C.
      • Li P.
      • Zhou G.
      • Chen S.
      • Tang Z.
      Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia.
      ). Furthermore, upregulated pyruvate dehydrogenase kinase isoform 1, which is a member of the HIF-1α target gene, acts to suppress pyruvate dehydrogenase, involved in inhibiting the initiation of the “TCA cycle.” The downregulation of all rate-limiting enzymes of the “TCA cycle” (citrate synthase, isocitrate dehydrogenase, dihydrolipoyl dehydrogenase, miR-193b/457a/15b-OGDH, miR-210/338-SDHb, SDHa, and miR-457a/214-MDH1) indicates inhibition of aerobic metabolism (Table 1 and Fig. 5) (
      • Cerretelli P.
      • Gelfi C.
      Energy metabolism in hypoxia: Reinterpreting some features of muscle physiology on molecular grounds.
      ).
      Surprisingly, both the rate-limiting enzymes (pyruvate kinase, miR-181a/-338-hexokinase 1 [HK1], HK2, and phosphofructokinase 1) of “glycolysis/gluconeogenesis” and l-lactate dehydrogenase were significantly downregulated (Table 1 and Fig. 5), contrary to the accepted metabolic model of hypoxic stress. Such a lower level for anaerobic glycolysis was found at both transcriptional and protein levels in P. vachelli muscle induced by hypoxia. However, these results were different from that of other fishes. For example, hypoxia increased glycolysis-related enzyme activities and l-lactate dehydrogenase in the P. vachelli liver and D. rerio muscle, respectively (
      • Chen K.
      • Cole R.B.
      • Rees B.B.
      Hypoxia-induced changes in the zebrafish (Danio rerio) skeletal muscle proteome.
      ). It may be that the muscles in our treatment undergo more severe hypoxia, resulting in a cellular metabolic disorder.

      Energy Metabolism

      Mitochondrial Dysfunction

      “Oxidative phosphorylation” enables the generation of ATP in the mitochondrial inner membrane and is the main source of energy in eukaryotes. It is generally accepted that hypoxia or ischemia inhibits mitochondrial oxidative phosphorylation and reduces the overall efficiency of energy metabolism (
      • Cerretelli P.
      • Gelfi C.
      Energy metabolism in hypoxia: Reinterpreting some features of muscle physiology on molecular grounds.
      ). We found that all 34 DEPs of “oxidative phosphorylation” were downregulated (Fig. 2), indicating that the coupling reaction of ADP and inorganic phosphate to ATP through the respiratory chain was inhibited to a certain extent. It also reflected that reduced locomotion powered by skeletal muscle might pose a minimal threat to survival under hypoxia when swimming activity will be minimized.
      Meanwhile, we found that downregulated expression of the regulator of mitochondrial biogenesis (peroxisome proliferator–activated receptor gamma coactivator-related protein 1; elongation factor G, mitochondrial; lon protease homolog 1, mitochondrial; protein SERAC1) and the downstream mitochondrial-related functional gene/proteins (cytochrome c oxidase 1/2/4 [COX1/2/4], Mt-cytochome b [Mt-Cyb], miR-457a/15b-UQCRC2 [cytochrome b-c1 complex subunit 1/2, mitochondrial], UQCRC1, CYC1 [cytochrome c1, heme protein, mitochondrial], Cyb561a3 [cytochrome b ascorbate-dependent protein 3], miR-457a/338/15b/PC-7236-ATP5A1 [ATP synthase subunit alpha/beta, mitochondrial], and ATP5B), which suggested that mitochondria may be dysfunctional (Table 1 and Fig. 5). Such mitochondrial damage induced by hypoxia may lead to the production of reactive oxygen species (ROS) and the decline of the membrane potential of mitochondria, causing a metabolic disorder and apoptosis (
      • Wang D.
      • Du Y.
      • Xu H.
      • Pan H.
      • Wang R.
      Paeonol protects mitochondrial injury and prevents pulmonary vascular remodeling in hypoxia.
      ).

      Catabolism and Anabolism

      P. vachelli has been shown to have an acute response to hypoxia, including activation of catabolic capacity for more energy, and reduction of biosynthetic capacity to reduce energy expenditure in the liver. Similar results are found in P. vachelli muscles, for example, upregulation of peroxisome proliferators–activated receptor alpha/beta is involved in lipid metabolism and lipid oxidation of skeletal muscle in low oxygen environments (
      • O'Brien K.A.
      • Horscroft J.A.
      • Devaux J.
      • Lindsay R.T.
      • Steel A.S.
      • Clark A.D.
      • Philp A.
      • Harridge S.D.R.
      • Murray A.J.
      PPARα-independent effects of nitrate supplementation on skeletal muscle metabolism in hypoxia.
      ) (Table 1). Also, hypoxia induces the activation of riboflavin catabolism through upregulation of RFK as the key enzyme, which improves carbohydrate, lipid, and amino acid metabolism in low oxygen environments (
      • Wang Y.P.
      • Wei J.Y.
      • Yang J.J.
      • Gao W.N.
      • Wu J.Q.
      • Guo C.J.
      Riboflavin supplementation improves energy metabolism in mice exposed to acute hypoxia.
      ). In this study, we also found ample evidence of significantly increased catabolism in P. vachelli muscles. The transport-related enriched pathways (e.g., SNARE interactions in vesicular transport and ABC transporters) and catabolic pathways (e.g., pyrimidine metabolism, central carbon metabolism in cancer) contained a large proportion of upregulated genes (Fig. 2, A and C). Furthermore, we found a significant decrease in intermediate metabolites in catabolism in metabolome KEGG analysis (e.g., metabolic pathways, glutathione, purine and linoleic acid metabolism), which indicated that the consumption of these intermediate metabolites under hypoxia exceeds their synthesis (Fig. 2C).
      Furthermore, the upregulated intermediate metabolites were a big part of anabolic pathways (e.g., biosynthesis of aminoacyl-tRNA, plant hormone, alkaloids derived from histidine, purine, ornithine, lysine, and nicotinic acid) in this study (Fig. 2C). This indicated that consumption of these intermediate metabolites was reduced and accumulated. In addition, we identified a large portion of downregulated genes/proteins enriched in amino acid biosynthesis pathways (e.g., aminoacyl-tRNA, ribosome, tryptophan, lysine, cysteine, and methionine). These results confirmed the downregulation of anabolism in fish muscle induced by hypoxia, especially in amino acid synthesis (
      • Hardy K.M.
      • Burnett K.G.
      • Burnett L.E.
      Effect of hypercapnic hypoxia and bacterial infection (Vibrio campbellii) on protein synthesis rates in the Pacific whiteleg shrimp, Litopenaeus vannamei.
      ) (Fig. 2, A and B).

      Muscle Function

      Our previous research found that in P. vachelli liver and brain tissues, the oxidative stress indices and antioxidant enzymatic activities were activated to varying degrees by hypoxia (
      • Zhang G.
      • Yin S.
      • Mao J.
      • Liang F.
      • Zhao C.
      • Li P.
      • Zhou G.
      • Chen S.
      • Tang Z.
      Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia.
      ). Similarly, we detected upregulated manganese superoxide dismutase and glutathione-S-transferase Mu 3/theta-3 mRNAs in P. vachelli muscles (Table 1). These might promote the inflammasome activation and trigger the ROS-induced apoptotic pathway (
      • Giraud-Billoud M.
      • Rivera-Ingraham G.A.
      • Moreira D.C.
      • Burmester T.
      • Castro-Vazquez A.
      • Carvajalino-Fernandez J.M.
      • Dafre A.
      • Niu C.
      • Tremblay N.
      • Paital B.
      • Rosa R.
      • Storey J.M.
      • Vega I.A.
      • Zhang W.
      • Yepiz-Plascencia G.
      • et al.
      Twenty years of the 'Preparation for Oxidative Stress' (POS) theory: Ecophysiological advantages and molecular strategies.
      ).

      Inflammatory Response

      Several of the KEGG pathways, enriched by transcriptomes with a large proportion of upregulated genes, were associated with human immune response (e.g., B-cell receptor signaling pathway and type II diabetes mellitus) and inflammation (e.g., nonalcoholic fatty liver disease, NOD-like receptor [NLR], and Fc epsilon RI signaling pathway). Of note, the NLR signaling pathway with 27 upregulated genes may offer insight into the molecular adaptations involved in hypoxia response (Fig. 2A). NLRs (NLR family CARD domain-containing protein 3/P3/P12, upregulated) detect the cytosolic presence of DAMP, mainly from mitochondria (e.g., ROS and mtDNA). Then NLRs drove the activation of NF-κB (interferon gamma receptor 1/IKBKG [NF-κB essential modulator]/NF-κB inhibitor beta/CHUK [inhibitor of nuclear factor kappa-B kinase subunit alpha], upregulated) and MAPK (MAP kinase-activated protein kinase 5/mitogen-activated protein kinase kinase kinase kinase 2/dual specificity mitogen-activated protein kinase kinase 4, upregulated), cytokine production, and apoptosis (
      • Chen Z.
      • Liu Y.
      • Sun B.
      • Li H.
      • Dong J.
      • Zhang L.
      • Wang L.
      • Wang P.
      • Zhao Y.
      • Chen C.
      Polyhydroxylated metallofullerenols stimulate IL-1beta secretion of macrophage through TLRs/MyD88/NF-kappaB pathway and NLRP(3) inflammasome activation.
      ,
      • Ko J.H.
      • Yoon S.O.
      • Lee H.J.
      • Oh J.Y.
      Rapamycin regulates macrophage activation by inhibiting NLRP3 inflammasome-p38 MAPK-NFkappaB pathways in autophagy- and p62-dependent manners.
      ). Alternatively, a different set of NLRs induced caspase-1 activation through the assembly of multiprotein complexes called inflammasomes. The activation of caspase-1 regulates the maturation of the proinflammatory cytokines interleukin 1β (IL-1β) and drives mitochondrial dysfunction and apoptosis (
      • Liu Q.
      • Zhang D.
      • Hu D.
      • Zhou X.
      • Zhou Y.
      The role of mitochondria in NLRP3 inflammasome activation.
      ,
      • Yi W.
      • Yiran L.
      • Bicheng L.
      Vicious circle of nlrp3 and mitochondrial damage plays central role in renal ischemia reperfusion injury.
      ) (Table 1).

      Apoptosis

      Studies have shown that the apoptosis of fish brain and cardiac cells in the hypoxic state is one of the main causes of “pond turnover” (
      • Xiao W.
      The hypoxia signaling pathway and hypoxic adaptation in fishes.
      ). Researchers have detected caspase-3/9 or apoptotic index and found that hypoxia stress can induce apoptosis of fish neural, liver, or cardiac cells, for example, Acipenser shrenckii (
      • Lu G.
      • Mak Y.T.
      • Wai S.M.
      • Kwong W.H.
      • Fang M.
      • James A.
      • Randall D.
      • Yew D.T.
      Hypoxia-induced differential apoptosis in the central nervous system of the sturgeon (Acipenser shrenckii).
      ), Gymnocephalus cernua, Platichthys flesus (
      • Tiedke J.
      • Thiel R.
      • Burmester T.
      Molecular response of estuarine fish to hypoxia: A comparative study with ruffe and flounder from field and laboratory.
      ), and M. amblycephala (
      • Sun S.
      • Xuan F.
      • Ge X.
      • Zhu J.
      • Zhang W.
      Dynamic mRNA and miRNA expression analysis in response to hypoxia and reoxygenation in the blunt snout bream (Megalobrama amblycephala).
      ). We first found that there was a significant gene involved in cell cycle arrest (cyclin-dependent kinase 10/12/13/cyclin I/14-3-3 protein zeta/delta, DNA replication licensing factor MCM4/7 downregulated, cyclin G2/death-inducer obliterator 1/mediator of DNA damage checkpoint protein 1/cyclin-dependent kinase inhibitor 1B/growth arrest and DNA damage–inducible protein GADD45 alpha upregulated) in P. vachelli muscles under hypoxia (Table 1), which may be due to the blocked energy metabolism and inflammatory response (
      • Skovira J.W.
      • Wu J.
      • Matyas J.J.
      • Kumar A.
      • Hanscom M.
      • Kabadi S.V.
      • Fang R.
      • Faden A.I.
      Cell cycle inhibition reduces inflammatory responses, neuronal loss, and cognitive deficits induced by hypobaria exposure following traumatic brain injury.
      ). C-fos (proto-oncogene c-Fos) and high mobility group protein B1/3 (HMGB1/3) were often the markers of cell apoptosis and play a significant role in hypoxia or ischemia-induced apoptosis, which were significantly upregulated in P. vachelli muscle (Fig. 6) (
      • Brunelle J.K.
      • Chandel N.S.
      Oxygen deprivation induced cell death: An update.
      ). “Apoptosis” was also significantly enriched in the KEGG pathway enrichment analysis of miRNA–mRNA pairs (Fig. 3A).
      The three pathways of apoptosis (extrinsic, mitochondrial, and endoplasmic reticulum stress) were enriched by abundant DEG/DEP (Fig. 6). For instance, in 145 negatively correlated miRNA–mRNA–protein pairs, HADHA (hydroxyacyl-CoA dehydrogenase subunit alpha) was downregulated and promoted cell apoptosis by decreasing the acylation of monolysocardiolipin into cardiolipin via raised miR-1-2 (
      • Taylor W.A.
      • Mejia E.M.
      • Mitchell R.W.
      • Choy P.C.
      • Sparagna G.C.
      • Hatch G.M.
      Human trifunctional protein alpha links cardiolipin remodeling to beta-oxidation.
      ). And that GRP78 is expressed in large quantities under low oxygen to maintain the stability of the endoplasmic reticulum and protect cells (
      • Reddy R.K.
      • Mao C.
      • Baumeister P.
      • Austin R.C.
      • Kaufman R.J.
      • Lee A.S.
      Endoplasmic reticulum chaperone protein GRP78 protects cells from apoptosis induced by topoisomerase inhibitors: Role of ATP binding site in suppression of caspase-7 activation.
      ). Upregulation of miR-457a led to endoplasmic reticulum stress and hypoxia-induced apoptosis by targeting GRP78 in P. vachelli muscles (Table 1 and Fig. 5).

      Muscle Dysfunction

      Several of the KEGG enriched by multiomics (p < 0.05) were associated with impaired contraction or dilation of the heart muscle (dilated cardiomyopathy and hypertrophic cardiomyopathy) and dyskinesia (Parkinson's and Huntington's disease), which suggested that hypoxia induces muscle dysfunction (Figs. 2 and 3A). In addition, “vascular smooth muscle contraction” and “calcium signaling pathway” contained abnormal muscular and calcium overload–related proteins that were significantly enriched by down DEP. For example, Atp2a1, a key regulator of muscle performance, contributed to calcium sequestration involved in muscular excitation/contraction (
      • Odermatt A.
      • Barton K.
      • Khanna V.K.
      • Mathieu J.
      • Escolar D.
      • Kuntzer T.
      • Karpati G.
      • MacLennan D.H.
      The mutation of Pro789 to Leu reduces the activity of the fast-twitch skeletal muscle sarco(endo)plasmic reticulum Ca2+ ATPase (SERCA1) and is associated with Brody disease.
      ). Krt8 helps to link the contractile apparatus to dystrophin at the sarcomere of striated muscle (
      • Ursitti J.A.
      • Lee P.C.
      • Resneck W.G.
      • McNally M.M.
      • Bowman A.L.
      • O'Neill A.
      • Stone M.R.
      • Bloch R.J.
      Cloning and characterization of cytokeratins 8 and 19 in adult rat striated muscle. Interaction with the dystrophin glycoprotein complex.
      ). In our results (Table 1 and Fig. 5), let-7b was upregulated, whereas its target genes/proteins Atp2a1 and Krt8 were significantly downregulated. This supported the idea that calcium overload and muscle dysfunction were linked to hypoxia adaptation in P. vachelli.
      To date, understanding of small RNAs in fishes is limited, and information regarding hypoxic miRNA markers are especially rare. Compared with other animals or cell models, 24 miRNAs reported previously as hypoxia responsive were identified (e.g., miR-15b, 133b, 143, 146b, 181a, 193b, 27b, 301b, 338, 152, 18a, 30a/d, 126, 199a, 210, 214, 29b, 25, let-7b) from 39 DEMIs in our results. These are thought to exert broad pleiotropic effects by targeting genes involved in cell cycle arrest, metabolism, apoptosis, cell survival, and mitochondrial function (
      • Bandara K.V.
      • Michael M.Z.
      • Gleadle J.M.
      MicroRNA biogenesis in hypoxia.
      ,
      • Guimbellot J.S.
      • Erickson S.W.
      • Mehta T.
      • Wen H.
      • Page G.P.
      • Sorscher E.J.
      • Hong J.S.
      Correlation of microRNA levels during hypoxia with predicted target mRNAs through genome-wide microarray analysis.
      ,
      • McCormick R.
      • Buffa F.M.
      • Ragoussis J.
      • Harris A.L.
      The role of hypoxia regulated microRNAs in cancer.
      ,
      • Pocock R.
      Invited review: Decoding the microRNA response to hypoxia.
      ). In addition, in fishes, the tendencies of miR-210, miR-193b, miR-181a expression among muscle were the same as the liver of P. vachelli (
      • Zhang G.
      • Yin S.
      • Mao J.
      • Liang F.
      • Zhao C.
      • Li P.
      • Zhou G.
      • Chen S.
      • Tang Z.
      Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia.
      ) or M. amblycephala (
      • Sun S.
      • Xuan F.
      • Ge X.
      • Zhu J.
      • Zhang W.
      Dynamic mRNA and miRNA expression analysis in response to hypoxia and reoxygenation in the blunt snout bream (Megalobrama amblycephala).
      ). For example, in negative miRNA–mRNA–protein pairs (supplemental Table S4), several key enzymes of the “TCA cycle,” for example, miR-193b-OGDH and miR-210-SDHb/SUCLA2, were detected in our conjoint analysis. Moreover, the same result as D. rerio that raised let-7b acts downstream of HIF-1α to repress cell proliferation through blocking cell cycle progression (
      • Huang C.X.
      • Chen N.
      • Wu X.J.
      • He Y.
      • Huang C.H.
      • Liu H.
      • Wang W.M.
      • Wang H.L.
      Zebrafish let-7b acts downstream of hypoxia-inducible factor-1alpha to assist in hypoxia-mediated cell proliferation and cell cycle regulation.
      ). Interestingly, several miRNA expressions in P. vachelli muscle were the opposite to that of liver (e.g., miR-27b, miR-143, miR-338), which could be due to different tissues or a metabolic disorder in muscle induced by hypoxia. For instance, data analysis of predicted multiomics pairs in silico (supplemental Table S4) suggested that the possible mechanism that GAPVD1, a protein required for efficient endocytosis (
      • Guillen R.X.
      • Beckley J.R.
      • Chen J.-S.
      • Gould K.L.
      CRISPR-mediated gene targeting of CK1δ/ε leads to enhanced understanding of their role in endocytosis via phosphoregulation of GAPVD1.
      ), promoted catabolism in muscles under hypoxia by downregulating miR-27b/143. Furthermore, upregulating miR-338 negatively targeted key enzymes of glycolysis (HK1) and “oxidative phosphorylation” (SDH, Atp5a1, and NDUFV1/S1) under hypoxic conditions in P. vachelli muscle (
      • Aschrafi A.
      • Kar A.N.
      • Natera-Naranjo O.
      • MacGibeny M.A.
      • Gioio A.E.
      • Kaplan B.B.
      MicroRNA-338 regulates the axonal expression of multiple nuclear-encoded mitochondrial mRNAs encoding subunits of the oxidative phosphorylation machinery.
      ). It is well known that many fish suppress locomotion during hypoxia as an energy-saving strategy. The resulting suppression or disorder of energy metabolism is probably beneficial.

      Conclusions

      We identified a large number of hypoxia-associated markers in the muscles of P. vachelli that were involved in diverse biological pathways such as HIF-1 signaling, energy metabolism, and muscle function. Compared with the P. vachelli livers, we found that muscles experienced significant hypoxia-associated changes, and this was evident in discrete tissue-specific patterns present. For example, P. vachelli livers have increased anaerobic glycolysis, heme synthesis, erythropoiesis, and inhibited apoptosis when exposed to under 4 h of hypoxia, whereas the opposite was true for muscles.
      To maintain normal physical activity, fishes react spontaneously to acute hypoxia, by activated catabolic pathways to generate more energy, decreased biosynthesis to reduce energy consumption, and shifted from aerobic to anaerobic metabolic contributions. According to the results of omics analysis, we infer that hypoxia induced muscle dysfunction through impairing mitochondrial function, activating inflammasome, and apoptosis (
      • Zhou R.
      • Yazdi A.S.
      • Menu P.
      • Tschopp J.
      A role for mitochondria in NLRP3 inflammasome activation.
      ). The hypoxia-induced mitochondrial dysfunction enhanced ROS generation and apoptosis, further triggering IL-1β production via the inflammasome activation (
      • Pomerantz B.J.
      • Reznikov L.L.
      • Harken A.H.
      • Dinarello C.A.
      Inhibition of caspase 1 reduces human myocardial ischemic dysfunction via inhibition of IL-18 and IL-1beta.
      ). In turn, IL-1β further impaired mitochondrial function and apoptosis by suppressing downstream mitochondrial biosynthesis–related proteins, resulting in a vicious circle between inflammasome activation and mitochondrial dysfunction (
      • Yi W.
      • Yiran L.
      • Bicheng L.
      Vicious circle of nlrp3 and mitochondrial damage plays central role in renal ischemia reperfusion injury.
      ) (Fig. 7). We hope that this study provides fundamental data on analytical methods as well as biological and mechanical insights into future research for environmental stress–induced hypoxia.
      Figure thumbnail gr7
      Fig. 7A schematic model of acute reactions in Pelteobagrus vachelli muscles caused by acute hypoxia. For abbreviations and explanations, please see the text.

      Data Availability

      All data generated or analyzed during this study are included in this published article and its supplemental data files. The raw data from the high-throughput sequencing data were deposited at the National Center for Biotechnology Information Sequence Read Archive (IDs: SUB7714154 and SUB7765255). The MS proteomics data have been deposited to the ProteomeXchange Consortium via the iProX partner repository with the dataset identifier PXD020425. Metabolome data related to this study have been publicly available on Metabolight (ID: MTBLS1888).

      Supplemental data

      This article contains supplemental data.

      Conflict of interest

      The authors declare no competing interests.

      Acknowledgments

      This study was supported by the National Natural Science Foundation of China (grant nos.: 32102754 and 32102760), National Key R & D Program of China (grant no: 2018YFD0900301), The Creation Project of Major New Species of Agriculture in Jiangsu Province (grant no: PZCZ201742), The “JBGS” Project of Seed Industry Revitalization in Jiangsu Province [grant no: JBGS(2021)034], Natural Science Foundation of the Higher Education Institutions of Jiangsu Province, China (grant nos.: 18KJB240001 and 20KJD240001), Qingchuang Science and Technology Support Program of Shandong Provincial College, and Natural Science Foundation of Shandong Province (ZR2019QC002). This work was supported in part by [email protected] administered by Innovation and Technology Commission.

      Author contributions

      S. Y. methodology; D. Y., Y. L., Y. Z., J. C., J. Z., K. Z., J. J., T. W., and Y. J. software; D. Y., Y. L., Y. Z., J. C., K. Z., J. J., T. W., and Y. J. validation; K. Z., J. J., T. W., and Y. J. formal analysis; D. Y., Y. L., Y. Z., J. C., K. Z., J. J., T. W., J. J., and Y. J. investigation; J. L. and G. Z. data curation; J. L. and G. Z. writing–original draft; D. Y, Y. L., Y. Z., J. C., K. Z., J. J., T. W., Y. J., and S. Y. writing–review & editing.

      Supplemental Data

      Figure thumbnail figs1
      Supplemental Figure S1Hierarchical clustering of DEG, DEMI, and DEP across six libraries and neg-DEM, and pos-DEM across 20 samples.

      References

        • Mandic M.
        • Ramon M.L.
        • Gracey A.Y.
        • Richards J.G.
        Divergent transcriptional patterns are related to differences in hypoxia tolerance between the intertidal and the subtidal sculpins.
        Mol. Ecol. 2015; 23: 6091-6103
        • Abdelrahman H.
        • ElHady M.
        • Alcivar-Warren A.
        Aquaculture genomics, genetics and breeding in the United States: Current status, challenges, and priorities for future research.
        BMC Genomics. 2017; 18: 191
        • Sappal R.
        • Fast M.
        • Purcell S.
        • MacDonald N.
        • Stevens D.
        • Kibenge F.
        • Siah A.
        • Kamunde C.
        Copper and hypoxia modulate transcriptional and mitochondrial functional-biochemical responses in warm acclimated rainbow trout (Oncorhynchus mykiss).
        Environ. Pollut. 2016; 211: 291-306
        • Zhong X.P.
        • Dan W.
        • Zhang Y.B.
        • Gui J.F.
        Identification and characterization of hypoxia-induced genes in Carassius auratus blastulae embryonic cells using suppression subtractive hybridization.
        Comp. Biochem. Physiol. B Biochem. Mol. Biol. 2009; 152: 161-170
        • Zhu C.D.
        • Wang Z.H.
        • Yan B.
        Strategies for hypoxia adaptation in fish species: A review.
        J. Comp. Physiol. B. 2013; 183: 1005-1013
        • Gracey A.Y.
        • Troll J.V.
        • Somero G.N.
        Hypoxia-induced gene expression profiling in the euryoxic fish Gillichthys mirabilis.
        Proc. Natl. Acad. Sci. U. S. A. 2001; 98: 1993-1998
        • Lardon I.
        • Eyckmans M.
        • Vu T.N.
        • Laukens K.
        • Boeck G.D.
        • Dommisse R.
        1H-NMR study of the metabolome of a moderately hypoxia-tolerant fish, the common carp (Cyprinus carpio).
        Metabolomics. 2013; 9: 1216-1227
        • Everett M.V.
        • Antal C.E.
        • Crawford D.L.
        The effect of short-term hypoxic exposure on metabolic gene expression.
        J. Exp. Zool. A Ecol. Genet. Physiol. 2012; 317: 9-23
        • Ju Z.
        • Wells M.C.
        • Heater S.J.
        • Walter R.B.
        Multiple tissue gene expression analyses in Japanese medaka (Oryzias latipes) exposed to hypoxia.
        Comp. Biochem. Physiol. C Toxicol. Pharmacol. 2007; 145: 134-144
        • Chen N.
        • Chen L.P.
        • Zhang J.
        • Chen C.
        • Wei X.L.
        • Gul Y.
        • Wang W.M.
        • Wang H.L.
        Molecular characterization and expression analysis of three hypoxia-inducible factor alpha subunits, HIF-1α/2α/3α of the hypoxia-sensitive freshwater species, Chinese sucker.
        Gene. 2012; 498: 81-90
        • Wang H.
        • Huang C.
        • Chen N.
        • Zhu K.
        • Chen B.
        • Wang W.
        • Wang H.
        Molecular characterization and mRNA expression of HIF-prolyl hydroxylase-2 (phd2) in hypoxia-sensing pathways from Megalobrama amblycephala.
        Comp. Biochem. Physiol. B Biochem. Mol. Biol. 2015; 186: 28-35
        • Beck B.H.
        • Fuller S.A.
        • Li C.
        • Green B.W.
        • Zhao H.
        • Rawles S.D.
        • Webster C.D.
        • Peatman E.
        Hepatic transcriptomic and metabolic responses of hybrid striped bass (Morone saxatilis×Morone chrysops) to acute and chronic hypoxic insult.
        Comp. Biochem. Physiol. Part D Genomics Proteomics. 2016; 18: 1-9
        • Gong D.
        • Xu L.
        • Li W.
        • Shang R.
        • Liu S.
        Comparative analysis of liver transcriptomes associated with hypoxia tolerance in the gynogenetic blunt snout bream.
        Aquaculture. 2020; 523: 735163
        • Yang Y.
        • Fu Q.
        • Wang X.
        • Liu Y.
        • Zeng Q.
        • Li Y.
        • Gao S.
        • Bao L.
        • Liu S.
        • Gao D.
        • Dunham R.
        • Liu Z.
        Comparative transcriptome analysis of the swimbladder reveals expression signatures in response to low oxygen stress in channel catfish, Ictalurus punctatus.
        Physiol. Genomics. 2018; 50: 636-647
        • Dowd W.W.
        • Renshaw G.M.
        • Cech Jr., J.J.
        • Kultz D.
        Compensatory proteome adjustments imply tissue-specific structural and metabolic reorganization following episodic hypoxia or anoxia in the epaulette shark (Hemiscyllium ocellatum).
        Physiol. Genomics. 2010; 42: 93-114
        • Smith R.W.
        • Cash P.
        • Ellefsen S.
        • Nilsson G.E.
        Proteomic changes in the crucian carp brain during exposure to anoxia.
        Proteomics. 2009; 9: 2217-2229
        • Wulff T.
        • Jokumsen A.
        • Hojrup P.
        • Jessen F.
        Time-dependent changes in protein expression in rainbow trout muscle following hypoxia.
        J. Proteomics. 2012; 75: 2342-2351
        • Sun S.
        • Guo Z.
        • Fu H.
        • Zhu J.
        • Ge X.
        Integrated metabolomic and transcriptomic analysis of brain energy metabolism in the male Oriental river prawn (Macrobrachium nipponense) in response to hypoxia and reoxygenation.
        Environ. Pollut. 2018; 243: 1154-1165
        • Zhang G.
        • Zhao C.
        • Wang Q.
        • Gu Y.
        • Li Z.
        • Tao P.
        • Chen J.
        • Yin S.
        Identification of HIF-1 signaling pathway in Pelteobagrus vachelli using RNA-seq: Effects of acute hypoxia and reoxygenation on oxygen sensors, respiratory metabolism, and hematology indices.
        J. Comp. Physiol. B. 2017; 187: 931-943
        • Sun S.
        • Xuan F.
        • Ge X.
        • Zhu J.
        • Zhang W.
        Dynamic mRNA and miRNA expression analysis in response to hypoxia and reoxygenation in the blunt snout bream (Megalobrama amblycephala).
        Sci. Rep. 2017; 7: 12846
        • Zhang G.
        • Yin S.
        • Mao J.
        • Liang F.
        • Zhao C.
        • Li P.
        • Zhou G.
        • Chen S.
        • Tang Z.
        Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia.
        Sci. Rep. 2016; 6: 22907
        • Zhang G.
        • Zhang J.
        • Wen X.
        • Zhao C.
        • Zhang H.
        • Li X.
        • Yin S.
        Comparative iTRAQ-based quantitative proteomic analysis of Pelteobagrus vachelli liver under acute hypoxia: Implications in metabolic responses.
        Proteomics. 2017; 17: 17-18
        • Zhang G.
        • Li J.
        • Zhang J.
        • Liang X.
        • Zhang X.
        • Wang T.
        • Yin S.
        Integrated analysis of transcriptomic, miRNA and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in heterosis.
        Mol. Cell. Proteomics. 2019; 18: 1437-1453
        • Patro R.
        • Duggal G.
        • Love M.I.
        • Irizarry R.A.
        • Kingsford C.
        Salmon provides fast and bias-aware quantification of transcript expression.
        Nat. Methods. 2017; 14: 417-419
        • Mortazavi A.
        • Williams B.A.
        • McCue K.
        • Schaeffer L.
        • Wold B.
        Mapping and quantifying mammalian transcriptomes by RNA-Seq.
        Nat. Methods. 2008; 5: 621-628
        • Robinson M.D.
        • McCarthy D.J.
        • Smyth G.K.
        edgeR: A bioconductor package for differential expression analysis of digital gene expression data.
        Bioinformatics. 2010; 26: 139-140
        • Li W.
        • Yang Y.
        • Liu Y.
        • Liu S.
        • Li X.
        • Wang Y.
        • Zhang Y.
        • Tang H.
        • Zhou R.
        • Li K.
        Integrated analysis of mRNA and miRNA expression profiles in livers of Yimeng black pigs with extreme phenotypes for backfat thickness.
        Oncotarget. 2017; 8: 114787-114800
        • Zhao C.
        • Zhang G.
        • Yin S.
        • Li Z.
        • Wang Q.
        • Chen S.
        • Zhou G.
        Integrated analysis of mRNA-seq and miRNA-seq reveals the potential roles of sex-biased miRNA-mRNA pairs in gonad tissue of dark sleeper (Odontobutis potamophila).
        BMC Genomics. 2017; 18: 613
        • Chen Z.
        • Wen B.
        • Wang Q.
        • Tong W.
        • Guo J.
        • Bai X.
        • Zhao J.
        • Sun Y.
        • Tang Q.
        • Lin Z.
        • Lin L.
        • Liu S.
        Quantitative proteomics reveals the temperature-dependent proteins encoded by a series of cluster genes in thermoanaerobacter tengcongensis.
        Mol. Cell. Proteomics. 2013; 12: 2266-2277
        • Cao M.
        • Li C.
        • Liu Y.
        • Cai K.
        • Zhou X.
        Assessing urinary metabolomics in giant pandas using chromatography/mass spectrometry: Pregnancy-related changes in the metabolome.
        Front. Endocrinol. (Lausanne). 2020; 11: 215
        • Wen X.
        • Hu Y.
        • Zhang X.
        • Wei X.
        • Yin S.
        Integrated application of multi-omics provides insights into cold stress responses in pufferfish Takifugu fasciatus.
        BMC Genomics. 2019; 20: 563
        • Geng X.
        • Feng J.
        • Liu S.
        • Wang Y.
        • Arias C.
        • Liu Z.
        Transcriptional regulation of hypoxia inducible factors alpha (HIF-alpha) and their inhibiting factor (FIH-1) of channel catfish (Ictalurus punctatus) under hypoxia.
        Comp. Biochem. Physiol. B Biochem. Mol. Biol. 2014; 169: 38-50
        • Baptista R.B.
        • Souza-Castro N.
        • Almeida-Val V.M.
        Acute hypoxia up-regulates HIF-1alpha and VEGF mRNA levels in Amazon hypoxia-tolerant Oscar (Astronotus ocellatus).
        Fish Physiol. Biochem. 2016; 42: 1307-1318
        • Chen K.
        • Cole R.B.
        • Rees B.B.
        Hypoxia-induced changes in the zebrafish (Danio rerio) skeletal muscle proteome.
        J. Proteomics. 2013; 78: 477-485
        • Cerretelli P.
        • Gelfi C.
        Energy metabolism in hypoxia: Reinterpreting some features of muscle physiology on molecular grounds.
        Eur. J. Appl. Physiol. 2011; 111: 421-432
        • Wang D.
        • Du Y.
        • Xu H.
        • Pan H.
        • Wang R.
        Paeonol protects mitochondrial injury and prevents pulmonary vascular remodeling in hypoxia.
        Respir. Physiol. Neurobiol. 2019; 268: 103252
        • O'Brien K.A.
        • Horscroft J.A.
        • Devaux J.
        • Lindsay R.T.
        • Steel A.S.
        • Clark A.D.
        • Philp A.
        • Harridge S.D.R.
        • Murray A.J.
        PPARα-independent effects of nitrate supplementation on skeletal muscle metabolism in hypoxia.
        Biochim. Biophys. Acta Mol. Basis Dis. 2019; 1865: 844-853
        • Wang Y.P.
        • Wei J.Y.
        • Yang J.J.
        • Gao W.N.
        • Wu J.Q.
        • Guo C.J.
        Riboflavin supplementation improves energy metabolism in mice exposed to acute hypoxia.
        Physiol. Res. 2014; 63: 341-350
        • Hardy K.M.
        • Burnett K.G.
        • Burnett L.E.
        Effect of hypercapnic hypoxia and bacterial infection (Vibrio campbellii) on protein synthesis rates in the Pacific whiteleg shrimp, Litopenaeus vannamei.
        Am. J. Physiol. Regul. Integr. Comp. Physiol. 2013; 305: R1356-R1366
        • Giraud-Billoud M.
        • Rivera-Ingraham G.A.
        • Moreira D.C.
        • Burmester T.
        • Castro-Vazquez A.
        • Carvajalino-Fernandez J.M.
        • Dafre A.
        • Niu C.
        • Tremblay N.
        • Paital B.
        • Rosa R.
        • Storey J.M.
        • Vega I.A.
        • Zhang W.
        • Yepiz-Plascencia G.
        • et al.
        Twenty years of the 'Preparation for Oxidative Stress' (POS) theory: Ecophysiological advantages and molecular strategies.
        Comp. Biochem. Physiol. A Mol. Integr. Physiol. 2019; 234: 36-49
        • Chen Z.
        • Liu Y.
        • Sun B.
        • Li H.
        • Dong J.
        • Zhang L.
        • Wang L.
        • Wang P.
        • Zhao Y.
        • Chen C.
        Polyhydroxylated metallofullerenols stimulate IL-1beta secretion of macrophage through TLRs/MyD88/NF-kappaB pathway and NLRP(3) inflammasome activation.
        Small. 2014; 10: 2362-2372
        • Ko J.H.
        • Yoon S.O.
        • Lee H.J.
        • Oh J.Y.
        Rapamycin regulates macrophage activation by inhibiting NLRP3 inflammasome-p38 MAPK-NFkappaB pathways in autophagy- and p62-dependent manners.
        Oncotarget. 2017; 8: 40817-40831
        • Liu Q.
        • Zhang D.
        • Hu D.
        • Zhou X.
        • Zhou Y.
        The role of mitochondria in NLRP3 inflammasome activation.
        Mol. Immunol. 2018; 103: 115-124
        • Yi W.
        • Yiran L.
        • Bicheng L.
        Vicious circle of nlrp3 and mitochondrial damage plays central role in renal ischemia reperfusion injury.
        Nephrol. Dial. Transplant. 2015; 30: iii3-iii4
        • Xiao W.
        The hypoxia signaling pathway and hypoxic adaptation in fishes.
        Sci. China Life Sci. 2015; 58: 148-155
        • Lu G.
        • Mak Y.T.
        • Wai S.M.
        • Kwong W.H.
        • Fang M.
        • James A.
        • Randall D.
        • Yew D.T.
        Hypoxia-induced differential apoptosis in the central nervous system of the sturgeon (Acipenser shrenckii).
        Microsc. Res. Tech. 2005; 68: 258-263
        • Tiedke J.
        • Thiel R.
        • Burmester T.
        Molecular response of estuarine fish to hypoxia: A comparative study with ruffe and flounder from field and laboratory.
        PLoS one. 2014; 9e90778
        • Skovira J.W.
        • Wu J.
        • Matyas J.J.
        • Kumar A.
        • Hanscom M.
        • Kabadi S.V.
        • Fang R.
        • Faden A.I.
        Cell cycle inhibition reduces inflammatory responses, neuronal loss, and cognitive deficits induced by hypobaria exposure following traumatic brain injury.
        J. Neuroinflammation. 2016; 13: 299
        • Brunelle J.K.
        • Chandel N.S.
        Oxygen deprivation induced cell death: An update.
        Apoptosis. 2002; 7: 475-482
        • Taylor W.A.
        • Mejia E.M.
        • Mitchell R.W.
        • Choy P.C.
        • Sparagna G.C.
        • Hatch G.M.
        Human trifunctional protein alpha links cardiolipin remodeling to beta-oxidation.
        PLoS one. 2012; 7e48628
        • Reddy R.K.
        • Mao C.
        • Baumeister P.
        • Austin R.C.
        • Kaufman R.J.
        • Lee A.S.
        Endoplasmic reticulum chaperone protein GRP78 protects cells from apoptosis induced by topoisomerase inhibitors: Role of ATP binding site in suppression of caspase-7 activation.
        J. Biol. Chem. 2003; 278: 20915-20924
        • Odermatt A.
        • Barton K.
        • Khanna V.K.
        • Mathieu J.
        • Escolar D.
        • Kuntzer T.
        • Karpati G.
        • MacLennan D.H.
        The mutation of Pro789 to Leu reduces the activity of the fast-twitch skeletal muscle sarco(endo)plasmic reticulum Ca2+ ATPase (SERCA1) and is associated with Brody disease.
        Hum. Genet. 2000; 106: 482-491
        • Ursitti J.A.
        • Lee P.C.
        • Resneck W.G.
        • McNally M.M.
        • Bowman A.L.
        • O'Neill A.
        • Stone M.R.
        • Bloch R.J.
        Cloning and characterization of cytokeratins 8 and 19 in adult rat striated muscle. Interaction with the dystrophin glycoprotein complex.
        J. Biol. Chem. 2004; 279: 41830-41838
        • Bandara K.V.
        • Michael M.Z.
        • Gleadle J.M.
        MicroRNA biogenesis in hypoxia.
        Microrna. 2017; 6: 80-96
        • Guimbellot J.S.
        • Erickson S.W.
        • Mehta T.
        • Wen H.
        • Page G.P.
        • Sorscher E.J.
        • Hong J.S.
        Correlation of microRNA levels during hypoxia with predicted target mRNAs through genome-wide microarray analysis.
        BMC Med. Genomics. 2009; 2: 15
        • McCormick R.
        • Buffa F.M.
        • Ragoussis J.
        • Harris A.L.
        The role of hypoxia regulated microRNAs in cancer.
        Curr. Top. Microbiol. Immunol. 2010; 345: 47-70
        • Pocock R.
        Invited review: Decoding the microRNA response to hypoxia.
        Pflugers Arch. 2011; 461: 307-315
        • Huang C.X.
        • Chen N.
        • Wu X.J.
        • He Y.
        • Huang C.H.
        • Liu H.
        • Wang W.M.
        • Wang H.L.
        Zebrafish let-7b acts downstream of hypoxia-inducible factor-1alpha to assist in hypoxia-mediated cell proliferation and cell cycle regulation.
        Life Sci. 2017; 171: 21-29
        • Guillen R.X.
        • Beckley J.R.
        • Chen J.-S.
        • Gould K.L.
        CRISPR-mediated gene targeting of CK1δ/ε leads to enhanced understanding of their role in endocytosis via phosphoregulation of GAPVD1.
        Sci. Rep. 2020; 10: 6797
        • Aschrafi A.
        • Kar A.N.
        • Natera-Naranjo O.
        • MacGibeny M.A.
        • Gioio A.E.
        • Kaplan B.B.
        MicroRNA-338 regulates the axonal expression of multiple nuclear-encoded mitochondrial mRNAs encoding subunits of the oxidative phosphorylation machinery.
        Cell Mol. Life Sci. 2012; 69: 4017-4027
        • Zhou R.
        • Yazdi A.S.
        • Menu P.
        • Tschopp J.
        A role for mitochondria in NLRP3 inflammasome activation.
        Nature. 2011; 469: 221-225
        • Pomerantz B.J.
        • Reznikov L.L.
        • Harken A.H.
        • Dinarello C.A.
        Inhibition of caspase 1 reduces human myocardial ischemic dysfunction via inhibition of IL-18 and IL-1beta.
        Proc. Natl. Acad. Sci. U. S. A. 2001; 98: 2871-2876