Hypoxia Affects HIF-1/LDH-A Signaling Pathway by Methylation Modification and Transcriptional Regulation in Japanese Flounder (Paralichthys olivaceus)

Simple Summary With global climate change and increased aquaculture production, fishes in natural waters or aquaculture systems are easily subjected to hypoxic stress. However, our understanding about their responsive mechanisms to hypoxia is still limited. Japanese flounder (Paralichthys olivaceus) is a widely cultivated marine economical flatfish, whose hypoxic responsive mechanisms are not fully researched. In this study, responses to hypoxia were investigated at blood physiological, biochemical, hormonal, and molecular levels. Responsive mechanisms of the HIF-1/LDH-A signaling pathway in epigenetic modification and transcriptional regulation were also researched. These results are important for enriching the theory of environmental responsive mechanisms and guiding aquaculture. Abstract Japanese flounder (Paralichthys olivaceus) responsive mechanisms to hypoxia are still not fully understood. Therefore, we performed an acute hypoxic treatment (dissolved oxygen at 2.07 ± 0.08 mg/L) on Japanese flounder. It was confirmed that the hypoxic stress affected the physiological phenotype through changes in blood physiology (RBC, HGB, WBC), biochemistry (LDH, ALP, ALT, GLU, TC, TG, ALB), and hormone (cortisol) indicators. Hypoxia inducible factor-1 (HIF-1), an essential oxygen homeostasis mediator in organisms consisting of an inducible HIF-1α and a constitutive HIF-1β, and its target gene LDH-A were deeply studied. Results showed that HIF-1α and LDH-A genes were co-expressed and significantly affected by hypoxic stress. The dual-luciferase reporter assay confirmed that transcription factor HIF-1 transcriptionally regulated the LDH-A gene, and its transcription binding sequence was GGACGTGA located at −2343~−2336. The DNA methylation status of HIF-1α and LDH-A genes were detected to understand the mechanism of environmental stress on genes. It was found that hypoxia affected the HIF-1α gene and LDH-A gene methylation levels. The study uncovered HIF-1/LDH-A signaling pathway responsive mechanisms of Japanese flounder to hypoxia in epigenetic modification and transcriptional regulation. Our study is significant to further the understanding of environmental responsive mechanisms as well as providing a reference for aquaculture.


Introduction
Hypoxia is a global problem that affects the function and biodiversity of the aquatic ecosystem, where the health and well-being of aquatic organisms are often jeopardized [1][2][3]. When exposed to a hypoxic environment, fishes often show changes in macro physiology and micro molecules, such as swimming behavior, electrophysiology, hormonal profiles, biochemistry, and hypoxia-responsive factors [4][5][6][7]. More specifically, hypoxia The anticoagulant blood was used to determine the quantity of red blood cells (RBC) and white blood cells (WBC), and the concentration (g/L) of hemoglobin (HGB). Specifically, within 24 h of sampling, these indicators were detected after being shaken and mixed well by referring to Ni et al., 2016 and strictly following the instructions of BC-1800 automatic blood cell analyzer (Mindray, Shenzhen, China) [52].

Biochemical Indicators in Serum
Non-anticoagulant blood was centrifuged (6000× g, 10 min, 25 • C) to obtain blood serum and then stored at −80 • C. Here, seven biochemical parameters related to metabolism were measured, which included the enzymatic activities of lactate dehydrogenase (LDH), alkaline phosphatase (ALP), alanine aminotransferase (ALT), and the concentration of glucose (GLU), total cholesterol (TC), triglyceride (TG), and albumin (ALB). According to Ni et al., 2016, we used a BS-180 automatic biochemical analyzer (Mindray, China) to collect these data and strictly followed the instructions of related kits (Mindray, China) [52].

Hormonal Indicator in Serum
The concentration of cortisol (COR) in blood serum was detected using a GC-911 gamma radioimmunoassay counter (ZONKIA, Hefei, China) and Iodine [ 125 I] cortisol radioimmunoassay kit (CNNC, Beijing, China) following the manufacturer's instructions.
To examine the evolutionary relationships of Japanese flounder with other species, we downloaded the HIF-1α and LDH-A amino acid sequences of some species from NCBI (https://www.ncbi.nlm.nih.gov/, accessed on 24 April 2022). Phylogenetic trees were constructed using Molecular Evolutionary Genetics Analysis software (MEGA 7.0, https://www.megasoftware.net/, accessed on 24 April 2022) with Neighbor-Joining (NJ) method, in which topological stability of the trees was evaluated under 1000 bootstrap replications.
The hybridization experiment of labeled RNA probe and target mRNA in muscle tissue slices was conducted according to the method of Kanda et al., 2011 and [53,54]. A fluorescence microscope (ECHO RVL-100-G, San Diego, CA , USA) was used to take digital images.

Quantitative Real-Time PCR and Expression Analysis
TRIzol reagent (Invitrogen, Waltham, MA, USA) was used to extract total RNA from skeletal muscle. RNA concentration was measured using a nucleic acid analyzer Biodropsis BD-1000 (OSTC, Beijing, China), and its integrity was tested using agarose gel electrophoresis. Afterwards, a PrimeScript™ RT reagent kit with gDNA Eraser (Takara, Kusatsu, Japan) was used to obtain cDNA.

Transcriptional Regulation Detection of HIF-1 on the LDH-A Gene
The method of dual-luciferase reporter assay was used to verify the predicted transcriptional regulation relationship of HIF-1 on the LDH-A gene.
Secondly, HIF-1α and HIF-1β expression plasmids (pc3.1~HIF-1α, pc3.1~HIF-1β) were constructed using single endonuclease digestion (BamHI (NEB, USA)) and homologous recombination. Specifically, the coding sequences (CDS) of HIF-1α and HIF-1β genes were amplified using high fidelity PCR with cDNA as template. The same sequences at both ends of the pcDNA3.1(+) plasmid digestion site were added to the 5 ends of forward primer or reverse primer (Supplementary Table S1), respectively, in the highfidelity PCR for homologous recombination. The pure target sequence, which was obtained using the FastPure ® Gel DNA Extraction Mini Kit (Vazyme, Nanjing, China), underwent homologous recombination with the linearized pcDNA3.1(+) plasmid to construct the expression plasmid. The constructed plasmid was transformed into DH5α competent cells (Vazyme, Nanjing, China) and transferred to LB (Luria-Bertani) solid medium. After the single colony was expanded in LB liquid medium and sequenced correctly (Sangon Biotech, Shanghai, China), we extracted the expression plasmid from expanded LB liquid medium using the EndoFree Mini Plasmid Kit II (TIANGEN, Beijing, China). Similarly, the reporter plasmid (pGL~LDH-A) was constructed using the above methods with the pGL3-Basic plasmid and the promoter sequence of the LDH-A gene. The slight difference was that the pGL3-Basic plasmid was digested using double endonuclease (SacI, HindIII (NEB, Ipswich, MA, USA)), and the promoter sequence of the LDH-A gene was obtained using high fidelity PCR under genomic DNA as the template.
Thirdly, plasmids were transfected into human embryonic kidney 293T (HEK293T) cells. In detail, we resuscitated (37 • C) HEK293T cells, which were cryopreserved in liquid nitrogen. The cells were cultured to 3rd or 4th generation, whose growth was normal and steady, in DMEM/High Glucose culture medium (Servicebio, Wuhan, China) with 10% fetal bovine serum (FBS) (absin, Shanghai, China) under 37 • C with 5% CO 2 . Subsequently, these cells were placed in 24-well plates (Corning, New York, NY, USA) at the concentration of 10 5 cells per well. Plasmids were transfected into HEK293T cells with Xfect™ Transfection Reagent Kit (Takara, Kusatsu, Japan) when the confluence of cells reached 50-70%. PRL-TK plasmid (Promega, Madison, WI, USA), could express Ranilla luciferase, and was used as the control to detect transfection efficiency. Here, triplicates were conducted. The medium was changed after transfection for 4 h.
For identifying the specific binding sites of transcription factor HIF-1, we performed fragment deletion and directed bases mutation in the LDH-A gene promoter, and constructed fragment deletion reporter plasmids (pGL~Lf2, pGL~Lf1, pGL~Lf0) and base mutation reporter plasmids (pGL~Lmp, pGL~Lm2), similar to the above methods. Especially, the mutated sequences in base mutation reporter plasmids were obtained using fusion PCR, whose primer sequences are also shown in Table S1 (Supplementary).

DNA Methylation Detection and Analysis
To thoroughly understand the mechanism of acute hypoxia stress affecting biomolecules, we detected the methylated modification status of HIF-1α and LDH-A genes using bisulfite modification and sequenced methods from the perspective of epigenetics. In detail, the Marine Animals DNA Kit (TIANGEN, Beijing, China) was used to extract genomic DNA from skeletal muscle tissue. DNA concentration was measured using a nucleic acid analyzer Biodropsis BD-1000 (OSTC, Beijing, China), and its integrity was detected using agarose gel electrophoresis. The DNA was stored at −20 • C.
Bisulfite modified DNA was obtained using a BisulFlash TM DNA Modification Kit (Epigentek, Farmingdale, NY, USA). Subsequently, a TaKaRa EpiTaq TM HS (for bisulfitetreated DNA) Kit (Takara, Kusatsu, Japan) was used to perform methylation-specific PCR (MS-PCR) with bisulfite modified DNA as the template. In the MS-PCR, CpG islands and primers (Supplementary Table S1) were designed using the online MethPrimer design software (http://www.urogene.org/methprimer/, accessed on 24 April 2022). The MS-PCR products were separated by agarose gel electrophoresis, and then the target sequences were purified using a Midi Purification Kit (TIANGEN, Beijing, China). Purified DNA sequences were connected with pEASY-T1 vector (TransGen, Beijing, China), and transferred into Trans1-T1 phage resistant chemically competent cells (TransGen, Beijing, China). Three biological repetitions were performed, and about 10 clones in each individual were sequenced to calculate the methylated modification level. To evaluate the bisulfite modified efficiency and ensure the credibility of experimental results, we calculated the percentage of converted cytosines (excluding cytosines in CpG dinucleotides). The formula is as follows: the converted percentage = [number of converted cytosines (excluding cytosines in CpG dinucleotides)]/[total number of cytosines (excluding cytosines in CpG dinucleotides)] × 100.

Statistical Analysis
Data were expressed as mean ± standard error (M ± SE). To determine significant differences among different groups, we conducted a one-way ANOVA accompanied with the Duncan's post hoc tests under the premises of normal distribution and variance homogeneity. p value was set at 0.05. It is worth mentioning that not every fish, but random individuals, were detected in the molecular experiment (q-PCR and DNA methylation detection) under less intra-group (technical repetition) differences in physiological indicators. In addition, linear regression and correlation analyses were performed to explore linear relationships between two variables (mRNA relative expression levels and methylation levels; HIF-1α and LDH-A mRNA relative expression levels). In the linear regression analysis, R was used to determine the strength of the correlation between two factors under the premise of a p value less than 0.05. The absolute value of R greater than 0.8, greater than 0.3 and less than 0.8, and less than 0.3 indicated that two factors were correlated strongly, weakly, and not, respectively. Here, SPSS 22.0 software (SPSS Co. Ltd., Chicago, IL, USA) and OriginPro 8.0 software (OriginLab, Northampton, MA, USA) were used for statistical analyses and graph construction.

Results
3.1. Effect of Acute Hypoxia on Physiology, Biochemistry, and Hormones 3.1.1. Blood Physiology (RBC, HGB, and WBC) RBC, HGB, and WBC had different fluctuating trends when hypoxia stress continued for 24 h. RBC was not significantly impacted by hypoxia treatment ( Figure 1A). However, HGB first underwent a significant decline (i.e., values fell by 15.58% at 1 h and were significantly lower at 3 and 6 h) and later returned to levels like the control (p < 0.05, Figure 1B). WBC in the 12 h group (125.67 ± 3.78 10 9 /L) was nearly 1.6 times that of the control group (p < 0.05, Figure 1C). With the prolongation of exposure to hypoxia, LDH activity increased, and reached a significantly higher level at 24 h than the control group (p < 0.05, Figure 2A). The ALP activity increased continuously for the first 3 h, was significantly higher in the 3 h treatment group than the control group (p < 0.05), and thereafter recovered ( Figure 2B). However, the ALT activity did not change significantly during the 24 h experimental timeline ( Figure 2C linear relationships between two variables (mRNA relative expression levels and methylation levels; HIF-1α and LDH-A mRNA relative expression levels). In the linear regression analysis, R was used to determine the strength of the correlation between two factors under the premise of a p value less than 0.05. The absolute value of R greater than 0.8, greater than 0.3 and less than 0.8, and less than 0.3 indicated that two factors were correlated strongly, weakly, and not, respectively. Here, SPSS 22.0 software (SPSS Co. Ltd., Chicago, IL, USA) and OriginPro 8.0 software (OriginLab, Northampton, Ma, USA) were used for statistical analyses and graph construction.

Effect of Acute Hypoxia on Physiology, Biochemistry, and Hormones
3.1.1. Blood Physiology (RBC, HGB, and WBC) RBC, HGB, and WBC had different fluctuating trends when hypoxia stress continued for 24 h. RBC was not significantly impacted by hypoxia treatment ( Figure 1A). However, HGB first underwent a significant decline (i.e., values fell by 15.58% at 1 h and were significantly lower at 3 and 6 h) and later returned to levels like the control (p < 0.05, Figure  1B). WBC in the 12 h group (125.67 ± 3.78 10 9 /L) was nearly 1.6 times that of the control group (p < 0.05, Figure 1C). With the prolongation of exposure to hypoxia, LDH activity increased, and reached

Concentration in Blood Serum
The concentration of GLU peaked at 1.87 ± 0.28 mmol/L in the 6 h hypoxic group (p < 0.05) and later recovered ( Figure 3A). TC concentration followed a fluctuating trend (ranging from 6.27 ± 3.32 to 12.71 ± 0.73 mmol/L) across the experimental timeline ( Figure  3B), while TG concentration significantly increased at longer hypoxia treatments (e.g., 6 and 24 h; p < 0.05, Figure 3C). There was no significant difference in ALB concentration among the 6 treatment groups ( Figure 3D).
(A) Figure 2. Effects of acute hypoxia on enzymatic activities of lactate dehydrogenase (LDH) (A), alkaline phosphatase (ALP) (B), and alanine aminotransferase (ALT) (C) in blood serum. Different letters indicate significant differences (p < 0.05). There was no significant difference in ALT activity among the six treatment groups.

Concentration in Blood Serum
The concentration of GLU peaked at 1.87 ± 0.28 mmol/L in the 6 h hypoxic group (p < 0.05) and later recovered ( Figure 3A). TC concentration followed a fluctuating trend (ranging from 6.27 ± 3.32 to 12.71 ± 0.73 mmol/L) across the experimental timeline ( Figure 3B), while TG concentration significantly increased at longer hypoxia treatments (e.g., 6 and 24 h; p < 0.05, Figure 3C). There was no significant difference in ALB concentration among the 6 treatment groups ( Figure 3D).

Concentration in Blood Serum
The concentration of GLU peaked at 1.87 ± 0.28 mmol/L in the 6 h hypoxic group (p < 0.05) and later recovered ( Figure 3A). TC concentration followed a fluctuating trend (ranging from 6.27 ± 3.32 to 12.71 ± 0.73 mmol/L) across the experimental timeline ( Figure  3B), while TG concentration significantly increased at longer hypoxia treatments (e.g., 6 and 24 h; p < 0.05, Figure 3C). There was no significant difference in ALB concentration among the 6 treatment groups ( Figure 3D).  There was no significant difference in ALB concentration among the six treatment groups.

COR Hormone
COR concentration increased and then decreased, where it was significantly higher in the 3 h group than the control and 24 h groups (p < 0.05, Figure 4).

Genetic Structure and Phylogenetic Analysis of HIF-1α and LDH-A
The results of gene sequence analyses are shown in Table S2 (Supplementary) and Figure S2 (Supplementary). The exon number of the HIF-1α gene (Gene ID: 109628061) There was no significant difference in ALB concentration among the six treatment groups.

COR Hormone
COR concentration increased and then decreased, where it was significantly higher in the 3 h group than the control and 24 h groups (p < 0.05, Figure 4).

Genetic Structure and Phylogenetic Analysis of HIF-1α and LDH-A
The results of gene sequence analyses are shown in Table S2 (Supplementary) and Figure S2

Genetic Structure and Phylogenetic Analysis of HIF-1α and LDH-A
The results of gene sequence analyses are shown in Table S2 (Supplementary) and Figure S2 (Supplementary). The exon number of the HIF-1α gene (Gene ID: 109628061) and the LDH-A gene (Gene ID: 109638975) was 15 and 7, respectively. The protein sequence information is shown in Table S3 and Figure S3 (Supplementary), in which five different domains (HLH, PAS, PAC, coiled coil, low complexity) were found in the HIF-1α protein and only one (low complexity) in the LDH-A protein. Furthermore, α helices, β sheets, and loop were all found in their three-dimensional structures ( Figure 5).
Biology 2022, 11, x FOR PEER REVIEW 13 of 29 and the LDH-A gene (Gene ID: 109638975) was 15 and 7, respectively. The protein sequence information is shown in Table S3 and Figure S3 (Supplementary), in which five different domains (HLH, PAS, PAC, coiled coil, low complexity) were found in the HIF-1α protein and only one (low complexity) in the LDH-A protein. Furthermore, α helices, β sheets, and loop were all found in their three-dimensional structures ( Figure 5). The HIF-1α phylogenetic tree ( Figure 6A) was grouped into two diverse clades (fish and other advanced vertebrates). The phylogenetic relationship of LDH-A was similar to that of HIF-1α except African clawed frog (Xenopus laevis) was closer to the fish clade (Figure 6B). The HIF-1α phylogenetic tree ( Figure 6A) was grouped into two diverse clades (fish and other advanced vertebrates). The phylogenetic relationship of LDH-A was similar to that of HIF-1α except African clawed frog (Xenopus laevis) was closer to the fish clade ( Figure 6B). and the LDH-A gene (Gene ID: 109638975) was 15 and 7, respectively. The protein se quence information is shown in Table S3 and Figure S3 (Supplementary), in which five different domains (HLH, PAS, PAC, coiled coil, low complexity) were found in the HIF 1α protein and only one (low complexity) in the LDH-A protein. Furthermore, α helices β sheets, and loop were all found in their three-dimensional structures ( Figure 5). The HIF-1α phylogenetic tree ( Figure 6A) was grouped into two diverse clades (fish and other advanced vertebrates). The phylogenetic relationship of LDH-A was similar to that of HIF-1α except African clawed frog (Xenopus laevis) was closer to the fish clade (Fig  ure 6B).

The Colocalization and Expression of HIF-1α and LDH-A
HIF-1α and LDH-A mRNA were co-located near the nucleus in muscle fiber cells by D-ISH (Figure 7).

The Colocalization and Expression of HIF-1α and LDH-A
HIF-1α and LDH-A mRNA were co-located near the nucleus in muscle fiber cells by D-ISH (Figure 7).
The mRNA relative expressions of HIF-1α increased with extending exposure time to hypoxia and then returned to the same level as the control group at 24 h ( Figure 8A). Throughout the overall hypoxia process, expression at 6 h was significantly higher than most other time groups (except 12 h) and was nearly 6.5-fold higher than the basal levels (p < 0.05). Moreover, there were no significant differences among the 0, 1, 3, and 24 h hypoxic treatment groups ( Figure 8A). LDH-A mRNA relative expressions were significantly higher in the 6 h group than the control group (p < 0.05, Figure 8B). The mRNA relative expressions of HIF-1α increased with extending exposure time to hypoxia and then returned to the same level as the control group at 24 h ( Figure 8A). Throughout the overall hypoxia process, expression at 6 h was significantly higher than most other time groups (except 12 h) and was nearly 6.5-fold higher than the basal levels (p < 0.05). Moreover, there were no significant differences among the 0, 1, 3, and 24 h hypoxic treatment groups ( Figure 8A). LDH-A mRNA relative expressions were significantly higher in the 6 h group than the control group (p < 0.05, Figure 8B).
Given the potential regulation of HIF-1 on the LDH-A gene, regression and correlation analyses were performed between their mRNA expressions. Here, we found that there was a positive correlation (R 2 = 0.589, R = 0.767, p < 0.05) between HIF-1α and LDH-A expressions ( Figure 8C).

HIF-1 Transcriptionally Regulates the LDH-A Gene
The dual-luciferase reporter assay result ( Figure 9A) showed that HIF-1 significantly upregulated the expression of the LDH-A gene (p < 0.05, 4th vs. 5th). In order to confirm HIF-1β being necessary for HIF-1, we also measured the relative luciferase activity values of the 6th group, in which expression plasmid pc3.1~HIF-1β was not added. The value of the 6th group was significantly lower than that of the 4th group (p < 0.05, b vs. c). However, interestingly, it was significantly higher than that of the 5th group (p < 0.05, b vs. a). In Figure 9B, the value in Lf2 was significantly higher than that in Lf1 (p < 0.05, b vs. a), which proved that the HIF-1 binding site was between −2443 and −2303. Combined with the results of JASPAR online software, we mutated the predicted binding bases (GGAC-GTGA, −2343~−2336). The values decreased significantly in Lmp (a vs. b) and Lm2 (a vs. b) compared with Lp and Lf2, respectively (p < 0.05, Figure 9C). Given the potential regulation of HIF-1 on the LDH-A gene, regression and correlation analyses were performed between their mRNA expressions. Here, we found that there was a positive correlation (R 2 = 0.589, R = 0.767, p < 0.05) between HIF-1α and LDH-A expressions ( Figure 8C).

HIF-1 Transcriptionally Regulates the LDH-A Gene
The dual-luciferase reporter assay result ( Figure 9A) showed that HIF-1 significantly upregulated the expression of the LDH-A gene (p < 0.05, 4th vs. 5th). In order to confirm HIF-1β being necessary for HIF-1, we also measured the relative luciferase activity values of the 6th group, in which expression plasmid pc3.1~HIF-1β was not added. The value of the 6th group was significantly lower than that of the 4th group (p < 0.05, b vs. c). However, interestingly, it was significantly higher than that of the 5th group (p < 0.05, b vs. a). In Figure 9B, the value in Lf2 was significantly higher than that in Lf1 (p < 0.05, b vs. a), which proved that the HIF-1 binding site was between −2443 and −2303. Combined with the results of JASPAR online software, we mutated the predicted binding bases (GGACGTGA, −2343~−2336). The values decreased significantly in Lmp (a vs. b) and Lm2 (a vs. b) compared with Lp and Lf2, respectively (p < 0.05, Figure 9C).

Methylation Modification and mRNA Expression of HIF-1α and LDH-A Genes
There are 18 CpG dinucleotides in the sequence of DNA methylation status measurement ( Figure 10A,B) of the HIF-1α gene (located at −243 to −10). Some putative transcription factor binding sites (TFBSs) are located near or at these CpG dinucleotides. We marked and showed these TFBSs and corresponding transcription factors (SOX10, Sox17, Sox6, VDR, ETS1, ZNF354C, GATA2, and GATA3) in Figure 10B. The sequence of DNA methylation status measurement ( Figure 10C) of the LDH-A gene (located at −753 to −511) contains 7 CpG dinucleotides ( Figure 10D). Similar to HIF-1α, TFBSs and transcription factors are also noted, which include SPl1, ZNF354C, ETS1, ARNT::HIF1A (HIF-1), Arnt, Mycn, TFE3, Creb312, USF1, and USF2. The pGL3 and pc3.1 represent circular pGL3-Basic plasmid and circular pcDNA3.1 (+) plasmid, respectively, which are not linked to the exogenous sequence. Reporter plasmid pGL~LDH-A was constructed using the LDH-A gene promoter sequence and fragmented pGL3-Basic plasmid digested by double endonuclease (SacI and HindIII). Expression plasmid pc3.1~HIF-1α was constructed using the HIF-1α CDS sequence and fragmented pcDNA3.1 (+) digested by single endonuclease (BamHI). Expression plasmid pc3.1~HIF-1β was constructed using the same method. PRL-TK is the control plasmid expressing Renilla luciferase (Rluc) to characterize the transfection efficiency. Arabic numerals of abscissa in Figure 9A indicate different treatment groups. (B) The fragmentation deletion result. The high triangle icons represent the putative HIF-1 binding site, which is predicted using JASPAR online software. Lp, Lf2, Lf1, and Lf0 represent reporter plasmids connecting different length promoter sequences, which are shown by four light gray horizontal lines. These reporter plasmids (pGL~Lf2, pGL~Lf1, and pGL~Lf0) were connected by different length LDH-A gene promoter sequences and fragmented pGL 3-Basic digested by double endonuclease (SacI and HindIII), respectively. Plasmids of pGL~Lp and pGL are the same as pGL~LDH-A and pGL3 in Figure 9A, respectively. (C) The bases mutation result. We mutated one putative transcription factor binding site (−2336) to obtain two mutant plasmids (pGL~Lmp, pGL~Lm2), which were deleted putative HIF-1 binding sequences (−2343~−2336: GGACGTGA). The mutant plasmids were constructed using fragmented pGL3-Basic digested using double endonuclease (SacI and HindIII) and the mutant sequence obtained using high fidelity fusion PCR. Different letters indicate significant differences (p < 0.05). In the process of HIF-1α gene methylation level detection, agarose gel electrophoresis showed that the size of MS-PCR products was the same as the target amplified product (234 bp). We randomly selected 12 sequencing results of the HIF-1α gene to calculate bisulfite modification efficiency. Only 17 cytosines (55 cytosines outside 18 CpG dinucleotides per 234 bp) were not converted to thymine, indicating that bisulfite modification efficiency was 97.42%, which was very efficient ( Figure S4A shows part of the sequencing results). The overall methylation level was 33.58 ± 2.30% across all time groups ( Figure  11A,B). DNA methylation level insignificantly changed when hypoxic stimulation oc- In the process of HIF-1α gene methylation level detection, agarose gel electrophoresis showed that the size of MS-PCR products was the same as the target amplified product (234 bp). We randomly selected 12 sequencing results of the HIF-1α gene to calculate bisulfite modification efficiency. Only 17 cytosines (55 cytosines outside 18 CpG dinucleotides per 234 bp) were not converted to thymine, indicating that bisulfite modification efficiency was 97.42%, which was very efficient ( Figure S4A shows part of the sequencing results). The overall methylation level was 33.58 ± 2.30% across all time groups (Figure 11A,B). DNA methylation level insignificantly changed when hypoxic stimulation occurred within 3 h, but then significantly declined in the 6 h group compared to that in the control and 1 h groups (p < 0.05, Figure 11B). Thereafter, at 12 and 24 h, DNA methylation level increased to comparable levels as the control group. In addition, methylation levels of single CpG dinucleotides first increased, decreased, and then recovered in most CpG dinucleotides (15/18) (Supplementary Figure S5).  Figure S5). Regression/correlation analyses showed that HIF-1α expression was negatively correlated with its methylation level (R 2 = 0.741, R = −0.861, p < 0.05) ( Figure 11C). Regression/correlation analyses was also performed between HIF-1α expression and methylation level of single CpG dinucleotides (Table S4 in   In the 243 bp sequence of the LDH-A gene, similar to the HIF-1α gene, we also randomly selected 12 sequencing results (containing 53 cytosines outside 7 CpG dinucleotides) which showed that only 11 cytosines were not converted to thymine, so that bisulfite modification efficiency was achieved to 98.27%; again, the results indicated high modification efficiency. Methylation level (9.55 ± 1.29%) was lower than the HIF-1α gene in all time groups (Figure 12A,B). The DNA methylation level first increased insignificantly at 1 h (p > 0.05), but it was significant at 3 h (p < 0.05). Thereafter, it significantly declined (at 6 and 12 h) and slightly increased again (at 24 h). In addition, single CpG dinucleotide methylation levels were similar to that of the HIF-1α gene (Supplementary Figure S6).
No correlation was found between LDH-A mRNA expression and the 243 bp overall methylation level (R 2 = 0.032). Furthermore, expression was similarly uncorrelated with 7 single CpG dinucleotides' methylation levels (Table S5 in Supplementary).

Response to Hypoxia Stress
HGB is a special protein that transports oxygen in RBC, which is even related to the immune system [58]. In this study, a significant decrease in HGB, but no significant increase in RBC was detected, indicating a decrease in oxygen delivery. The significant rebound of HGB and the slight retreat of RBC meant that fish were gradually adapting themselves to the hypoxic environment. WBC takes part in the process of cell-mediated immune response and phagocytosis directly in fish [59]. The fluctuating variations represented immune responses being activated in this hypoxic environment. Similarly, hypoxia also caused physiological response or adaption in Wuchang Bream (Megalobrama amblycephala), largemouth bass (Micropterus salmoides), and cavefish (Astyanax mexicanus) [60][61][62]. In short, hypoxia generally induced blood physiological responses in fish.
LDH, a carbohydrate metabolism enzyme, is a biomarker of stress and tissue damage [63,64]. Its enzymatic activity reflects the ability of anaerobic metabolism in organisms. Regression/correlation analyses showed that HIF-1α expression was negatively correlated with its methylation level (R 2 = 0.741, R = −0.861, p < 0.05) ( Figure 11C). Regression/correlation analyses was also performed between HIF-1α expression and methylation level of single CpG dinucleotides (Table S4 in Figure 11D).
In the 243 bp sequence of the LDH-A gene, similar to the HIF-1α gene, we also randomly selected 12 sequencing results (containing 53 cytosines outside 7 CpG dinucleotides) which showed that only 11 cytosines were not converted to thymine, so that bisulfite modification efficiency was achieved to 98.27%; again, the results indicated high modification efficiency. Methylation level (9.55 ± 1.29%) was lower than the HIF-1α gene in all time groups ( Figure 12A,B). The DNA methylation level first increased insignificantly at 1 h (p > 0.05), but it was significant at 3 h (p < 0.05). Thereafter, it significantly declined (at 6 and 12 h) and slightly increased again (at 24 h). In addition, single CpG dinucleotide methylation levels were similar to that of the HIF-1α gene (Supplementary Figure S6).
No correlation was found between LDH-A mRNA expression and the 243 bp overall methylation level (R 2 = 0.032). Furthermore, expression was similarly uncorrelated with 7 single CpG dinucleotides' methylation levels (Table S5 in Supplementary).

Response to Hypoxia Stress
HGB is a special protein that transports oxygen in RBC, which is even related to the immune system [58]. In this study, a significant decrease in HGB, but no significant increase in RBC was detected, indicating a decrease in oxygen delivery. The significant rebound of HGB and the slight retreat of RBC meant that fish were gradually adapting themselves to the hypoxic environment. WBC takes part in the process of cell-mediated immune response and phagocytosis directly in fish [59]. The fluctuating variations represented immune responses being activated in this hypoxic environment. Similarly, hypoxia also caused physiological response or adaption in Wuchang Bream (Megalobrama amblycephala), largemouth bass (Micropterus salmoides), and cavefish (Astyanax mexicanus) [60][61][62]. In short, hypoxia generally induced blood physiological responses in fish.
LDH, a carbohydrate metabolism enzyme, is a biomarker of stress and tissue damage [63,64]. Its enzymatic activity reflects the ability of anaerobic metabolism in organisms. Increased LDH activities in blood serum indicated that anaerobic metabolism enhanced, which was an adaptation strategy to hypoxia. The effect of changing dissolved oxygen concentrations on LDH activities were also found in cichlid (Astronotus Ocellatus) and rat [65,66]. All these findings and the result of LDH-A expression increase in skeletal muscle suggested that LDH was affected by hypoxia stress. ALP and ALT are liver marker enzymes [67]. Liver disease has also been further proved to be related to hypoxia stress in mice [68,69]. As a result, the fluctuation of ALP and ALT activities implied that these fishes suffered from liver damage under hypoxia.
Increased GLU and altered electrolyte homeostasis in blood are regarded as organism responses to environmental stress [70]. GLU was metabolized to provide energy by the glycolysis process in fish [71]. In our study, the result of GLU concentration gradually increasing when exposed to hypoxia indicates that metabolic stress responses occurred in these fishes. COR is one of the principal stress hormones in teleost fishes, and plays a significant role in osmoregulation, growth, and reproduction [72,73]. The significant increase in COR meant that our fishes were indeed affected by acute hypoxia stress. Furthermore, higher variability at 24 h was likely related to the unstable state of fishes. The same as Japanese flounder in our study, COR was also affected by hypoxia and other environmental pressure factors in rainbow trout (Oncorhynchus mykiss) and brown trout (Salmo trutta) [74,75]. Nevertheless, it has been noted that the importance of GLU is minor according to the absence of a unique mineralocorticoid in fishes, and it is the COR that has both glucocorticoid and mineralocorticoid roles [73]. This could explain why GLU concentration was not significantly changed among all time groups except for the 6 h group. Both TC and TG are related to lipid metabolism. Their fluctuations meant that lipid metabolisms were unstable. Lipid metabolites being significantly influenced by hypoxia was also found in Selenka (Apostichopus japonicus) and mice [2,76]. These findings implied that hypoxia stress may generally affect the homeostasis of the synthesis and degradation of lipids and their derivatives in many species.
The changes in physiological, biochemical, and hormonal indicators showed that acute hypoxic stress caused the stressful responses of Japanese flounder.
Correlations between gene expression and methylation have been found in many species. For example, Goyal and Goyal (2019) found that hypomethylation was associated with upregulated expression [83]. A negative correlation between mRNA expression and DNA methylation was also discovered by Li et al., 2017 [84]. Similar negative correlations between gene expressions and methylated levels were also found in the cyp19a1a gene of Japanese flounder and the SNCA gene of humans (Homo sapiens) [85,86]. Overall, the above findings and our result of negative correlations between HIF-1α gene expression and methylation level further indicated the potential regulation of DNA methylation on gene expression. Attwood et al., 2002 concluded that DNA methylation did affect gene expression by preventing transcriptional activation [25]. However, it is puzzling that the expression of LDH-A was not correlated with its methylation level. We speculated that the 243 bp sequence, which we detected, was not the important regulatory related sequence. Thomson et al., 2022 found that ocular growth rates were related to single site methylation levels [87]. Therefore, combining with the strong negative correlations of -188 and -78 sites in the HIF-1α gene, we are speculating that the methylation status of key nucleotide sites may be more important than the overall methylation level of CpG islands.

Transcriptional Regulation of HIF-1 and Its Target Genes
The HREs in HIF-1 transcription factor target genes are relatively conserved like HIF-1α and LDH-A [88]. The phenomenon of HIF-1 regulating LDH-A gene expression was found in mice [19]. Kuo et al., 2017 discovered that HIF-1α upregulated LDH-A expression in human (Homo sapiens) hypoxic neuroblastoma cells [89]. Xie et al., 2009 also pointed out that the LDH-A gene was one of the HIF-1α target genes in human cells [90]. Both these reports and our results of D-ISH, i.e., the positive correlation between HIF-1α and LDH-A gene expression, all implied that HIF-1 may transcriptionally regulate LDH-A gene expression in Japanese flounder. However, the importance of HIF-1β in HIF-1 and the specific binding sites of HIF-1 on the LDH-A gene still need to be further explored in Japanese flounder.
The dual-luciferase reporter assay proved the above conjecture ( Figure 9A). In addition, we also confirmed that HIF-1β was a synergistic factor for HIF-1α transcriptional regulation. The significant difference between the 5th and 6th treatments in Figure 9A was regarded as the endogenous HIF-1β interference in HEK293T cells. In brief, the results implied that the LDH-A gene was the HIF-1 target gene in Japanese flounder. Furthermore, the binding site of the HIF-1 transcription factor was located at −2343~−2336 (GGACGTGA) according to the results of fragment deletion and directed bases mutation experiment. Therefore, the sequence in sites of −2343~−2336 (GGACGTGA) may be the most important sequence for LDH-A gene sensing hypoxia. Additionally, its mutation will possibly affect transcriptional activation, and then affect a series of downstream cascade amplification effects of the signal pathway and is unfavorable to organisms or even affect the normal life activities. Here, we also conclude that not all predicted HRE sequences are the binding sites of HIF-1 transcription factors, which need to be identified by experiments rather than prediction. In conclusion, transcriptional regulation of HIF-1 on LDH-A may be common in organisms, but its specific nucleotide binding sites in many species need to be further studied.
Hypoxia exposure can result in large changes in gene expression, which is considered to be a response and potentially enhances hypoxic survival [91]. It has been widely reported that the target gene of HIF-1 (consisting of inducible HIF-1α and constitutive HIF-1β) includes GLUT-1, Epo, VEGF, and Hsp27 [12,60,92]. The results of HIF-1α being affected by hypoxia were found in largemouth bass (Micropterus salmoides) [60], oscar (Astronotus ocellatus) [93], Atlantic croaker (Micropogonias undulatus) [94], and Wuchang Bream (Megalobrama amblycephala) [62]. The changes in HIF-1α could influence the transcription of its target genes; thus, affect the whole metabolic network. Specifically, HIF-1, as a transcription factor, can recognize and bind specific promoters of target genes, recruit RNA polymerase, turn on transcription, and then improve target gene transcription levels, resulting in a series of downstream events. Luo and Wang (2018) also concluded that HIF was related to many human diseases in physiology and pathogenesis [95]. Therefore, the important regulatory role of HIF-1 response to hypoxic stress in many species needs to be further studied in the future.
The signal networks composed of hypoxia-related signaling pathways were activated by hypoxia to adapt themselves to the hypoxic environment in cells and organisms. These pathways were mainly related to respiration, metabolism, and survival [96]. Therefore, the hypoxia signaling networks in Japanese flounder need to be further studied in the future.

Conclusions
Hypoxia caused stress responses in physiology, biochemistry, hormone, mRNA expression, and DNA methylation modification of Japanese flounder. We studied the epigenetic modification and transcriptional regulation responsive mechanisms of HIF-1α (a typical inducible hypoxia-responsive factor) and its target gene LDH-A. In the study, it was concluded that hypoxia affected HIF-1 by methylated modification, and HIF-1 could transcriptionally regulate its target gene LDH-A by binding to a hypoxia response element (HRE) of GGACGTGA located at −2343~−2336, and then LDH-A affected the downstream cascade amplification, and further influenced the physiological phenotype of Japanese flounder. These physiological and molecular changes indicated that environmental factors caused various coordinated responses at multiple levels in organisms. These findings will provide a reference for the role of the environment on biomolecules and biomolecular interactions.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology11081233/s1, Table S1: Primer sequences for methylationspecific PCR (MS-PCR), quantitative real-time PCR (q-PCR), double in situ hybridization (D-ISH) and dual-luciferase reporter assay; Table S2: The gene sequence analysis results of HIF-1α and LDH-A; Table S3: The protein sequence information of HIF-1α and LDH-A; Table S4: Linear regression analysis results of expression and 18 single CpG dinucleotides methylation level in HIF-1α gene; Table S5: Linear regression analysis results of expression and 7 single CpG dinucleotides methylation level in LDH-A gene; Figure S1: The melt curves of primers for three genes HIF-1α (A), LDH-A (B), 18S (C) in q-PCR; Figure S2: The gene sequence analysis results of HIF-1α and LDH-A; Figure S3: Domain analyses of HIF-1α and LDH-A protein; Figure S4: The part DNA sequences of HIF-1α gene (A) and LDH-A gene (B) modified by bisulfite; Figure S5: The single CpG dinucleotides methylation level of HIF-1α gene in the six treatment groups of Japanese flounder (Paralichthys olivaceus); Figure  S6: The single CpG dinucleotides methylation level of LDH-A gene in the six hypoxia treatment groups of Japanese flounder. Data Availability Statement: Data used to support the findings from this study are available from the corresponding author upon request.