Front. Mar. Sci. Frontiers in Marine Science Front. Mar. Sci. 2296-7745 Frontiers Media S.A. 10.3389/fmars.2021.736362 Marine Science Original Research miRNA–mRNA Integrative Analysis Reveals the Roles of miRNAs in Hypoxia-Altered Embryonic Development- and Sex Determination-Related Genes of Medaka Fish Lai Keng Po 1 2 3 * Tam Nathan Yi Kan 2 Chen Yuelong 4 Leung Chi Tim 2 Lin Xiao 5 Tsang Chau Fong 2 Kwok Yin Cheung 2 Tse William Ka Fai 6 Cheng Shuk Han 7 Chan Ting Fung 4 Kong Richard Yuen Chong 2 3 * 1Laboratory of Environmental Pollution and Integrative Omics, Guilin Medical University, Guilin, China 2Department of Chemistry, City University of Hong Kong, Kowloon, Hong Kong SAR, China 3State Key Laboratory of Marine Pollution, City University of Hong Kong, Kowloon, Hong Kong SAR, China 4State Key Laboratory of Agrobiotechnology, School of Life Sciences, The Chinese University of Hong Kong, Shatin, Hong Kong SAR, China 5Department of Psychiatry, Icahn School of Medicine at Mount Sinai, New York, NY, United States 6Center for Promotion of International Education and Research, Faculty of Agriculture, Kyushu University, Fukuoka, Japan 7Department of Biomedical Sciences, City University of Hong Kong, Kowloon, Hong Kong SAR, China

Edited by: Youji Wang, Shanghai Ocean University, China

Reviewed by: Jae-Sung Rhee, Incheon National University, South Korea; Bindhu Paul, Amrita Vishwa Vidyapeetham, India

 *Correspondence: Keng Po Lai, kengplai@cityu.edu.hk Richard Yuen Chong Kong, bhrkong@cityu.edu.hk

This article was submitted to Aquatic Physiology, a section of the journal Frontiers in Marine Science

 21 01 2022 2021 8 736362 08 07 2021 14 12 2021 Copyright © 2022 Lai, Tam, Chen, Leung, Lin, Tsang, Kwok, Tse, Cheng, Chan and Kong. 2022 Lai, Tam, Chen, Leung, Lin, Tsang, Kwok, Tse, Cheng, Chan and Kong

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Recent studies have shown hypoxia to be an endocrine disruptor that impairs sex differentiation and reproductive function, leading to male-biased F1 populations in fish. However, the molecular mechanisms through which hypoxia alters fish sex differentiation and therefore sex ratios remain poorly understood. In order to understand the potential role of miRNAs in mediating hypoxia-altered sex determination and differentiation in fish, we conducted small RNA sequencing and transcriptome sequencing on marine medaka (Oryzias melastigma) embryos that were exposed to hypoxia (2.0 ± 0.2 mg O2 L–1) for 40 h (encompassing a critical window of sex determination). We identified dysregulated miRNAs and mRNAs in the hypoxia-exposed embryo, and bioinformatic analysis of the integrative small RNA sequencing and transcriptome sequencing results revealed hypoxia to cause alterations of genes related to embryonic development through miRNA regulation. Importantly, we have identified miRNA-mRNA pairs that were reported to play roles in gonad development (novel miR-145-col9a3 and novel miRNA-94- arid5b), in sex hormone response (novel miRNA-210-ca2, novel miRNA-106-nr2f2, nbr-miR-29c-nr4a1, and ola-miR-92b-akr1d1), and in sex characteristic development (novel miRNA-145-mns1, nle-miR-20-sord, and ipu-miR-219b-abcc8). Our findings highlighted the possible roles of miRNA–mRNA in regulation of embryonic development and sex determination in response to hypoxic stress.

 hypoxia fish sequencing miRNA bioinformatics Introduction

Hypoxia is a widespread and pressing environmental concern in aquatic habitats, causing severe habitat damage and major disruption to aquatic ecosystems worldwide. The occurrence of hypoxia in water bodies has increased globally over the past 30 years due to escalating eutrophication and organic pollution. Presently, over 500 hypoxic areas (<2 mg O2 L–1) spanning hundreds of thousands of square kilometers have been reported worldwide (Thrash et al., 2017), and is likely to worsen in the future, as a result of global warming and rapid coastal development, thereby threatening the sustainability of natural populations (Pörtner and Peck, 2010). Compared with mammals, fishes possess an exceptional range of sex determination and differentiation processes, which are modulated by a variety of environmental factors such as population density, social behaviors, pH, and dissolved oxygen (Baroiller et al., 2009; Reddon and Hurd, 2013; Yamamoto et al., 2014).

Previous studies from our group have shown that hypoxia is an endocrine disruptor that impairs reproductive activities and affects sexual differentiation in fish, leading to male-biased F1 generations (Shang et al., 2006; Cheung et al., 2014). Similar findings on sex ratio alterations have been observed in other studies which showed that the forkhead transcription factor (foxl2) and doublesex and mab3-related transcription factor 1 (dmrt1) were downregulated following exposure to hypoxia (Pelley, 2006; Baroiller et al., 2009). Importantly, a largescale field study from the Gulf of Mexico, one of the largest hypoxic dead zones in the world, reported similar endocrine disruption phenomena − extensive reproductive disruption, ovarian masculinization and male-biased Atlantic croaker sex ratios − occurring in natural fish populations (Thomas and Rahman, 2012). DNA methylation of the gonadal cyp19a (aromatase) gene has been linked to temperature-dependent sex ratio shifts in European sea bass (Navarro-Martín et al., 2011), suggesting that epigenetic modifications in response to environmental changes may also play a role in sex determination. Additionally, several recent studies have reported regulatory roles of miRNAs in sex differentiation. For example, let-7 and miR-21 have been shown to regulate egg development in rainbow trout (Ma et al., 2012), and miR-141 and miR-429 have been found to play crucial roles in regulating testis development and spermatogenesis in yellow catfish (Jing et al., 2014). Distinct subsets of miRNAs have also been observed in both male and female Nile tilapia embryos (Eshel et al., 2014) and gonads (Wang et al., 2016), indicating potential regulatory roles for miRNAs in the sexual differentiation of fish. Overall, the abovementioned studies strongly implicate miRNAs as a new and important class of effector molecules controlling gonadal development and sexual differentiation in vertebrates.

Marine medaka Oryzias melastigma (synonyms: Oryzias dancena) is a model marine fish species for ecotoxicological studies (Lau et al., 2014). The sex-determining gene DMY is absent in the O. melastigma (Juanchich et al., 2013). Genetic studies using a progeny test of sex-reversed individuals indicated the presence of XX/XY sex-determination system absent in the O. melastigma (Bouchareb et al., 2017). In addition, SRY-related HMG-Box gene 3 (sox3) has been co-opted as the male sex determination gene in the O. melastigma (Gay et al., 2018). Our group has previously identified several gonad-specific miRNAs in marine medaka, O. melastigma, as well as specific miRNAs from medaka embryos that are responsive to hypoxia (Feng et al., 2014). Likewise, differentially expressed miRNAs was found to play important roles during fish ovarian development (Kang et al., 2004). A genome-wide study identified novel ovarian-predominant miRNAs in the medaka fish (Mishima et al., 2008). And, miR-202 was reported to control female fecundity by regulating medaka oogenesis (Fagegaltier et al., 2014), suggesting that tissue-specific miRNAs may play specific roles in gonadal functions. Using marine medaka embryos as a model, together with integrative omics analysis, including small-RNA sequencing and transcriptome sequencing, we investigated the role of specific hypoxia-responsive miRNAs in directly targeting sex determination genes. Our results identify miRNA–mRNA pairs that could be important for controlling sexual differentiation, resulting in hypoxia-altered sex ratios in fish.

 Materials and Methods Fish Embryo Maintenance and Hypoxic Exposure

All animal research procedures were approved by the Animal Ethics Committee on Research Experiments involving Animal Subjects (A-0244) of the City University of Hong Kong. Oryzias melastigma embryos used in our experiments were obtained from State Key Laboratory of Marine Pollution, City University of Hong Kong. Embryos within 1 day old (< 24 h after fertilization) were collected from parental marine medakas that were maintained in seawater with salinity of 3.5% under optimal growth and breeding conditions (5.8 mg O2 L–1, 28 ± 2°C, pH 7.2 in a 14-h light: 10-h dark cycle). The collected embryos were randomly selected under a dissecting microscope. Filaments attached to the clustered embryos were carefully removed via constant rolling of the embryos on a net with fingertips and pulling of forceps until the embryo was fully separated and clean. Abnormal or unfertilized embryos were removed from the batch. The selected embryos were then distributed into 12 separate 50-mL beakers, each holding 30 embryos. Six beakers of embryos were exposed to normoxia (5.8 ± 0.2 mg O2 L–1; outside the hypoxic chamber) and six beakers of embryos were exposed to hypoxia (2.0 ± 0.2 mg O2 L–1; inside the hypoxic chamber) continuously for 40 h. The desired dissolved oxygen (DO) level was achieved by a constant flow of premixed air and nitrogen (0.5% oxygen) inside a hypoxic chamber. The DO level was measured and monitored using a DO meter (YSI model 580). After the exposure period, the DO level of the hypoxic water was measured again to ensure DO levels were maintained at 2.0 ± 0.2 mg O2 L–1. Water temperature was maintained at 26 ± 2°C by placing the hypoxic chamber inside an incubator.

 Embryo Survivability and Hatchability

After the exposure to normoxic or hypoxic conditions, embryos from each condition were transferred back to normoxic condition. The number of hatched eggs and dead eggs was recorded daily for 4 weeks. Temperature was maintained at 26 ± 2°C under a 14:10 light:dark photoperiod. Statistical analyses were performed using the paired t-test in GraphPad Prism 9.1.0 (GraphPad Software).

 RNA Sample Preparation and Small RNA Sequencing

After the normoxic or hypoxic exposure, one beaker of embryos from each condition was used as one biological replicate for the RNA preparation, and there were three replicates per condition. Total RNA was extracted immediately using the mirVanaTM miRNA isolation kit (Applied Biosystems) following the manufacturer’s instructions. RNA quality was assessed using the Agilent 2100 Bioanalyzer system and samples with an RNA Integrity Number (RIN) greater than 8 were used to construct the small RNA and complementary DNA (cDNA) libraries (Supplementary Table 1). Five μg of total RNA were used from each sample. Short RNA transcripts (18–30 nucleotides long) were resolved and isolated using polyacrylamide gel electrophoresis (PAGE) gels. The isolated small RNA molecules were ligated with 3′ and 5′ adapters, and single-stranded cDNA was synthesized using SuperScript II Reverse Transcriptase (Invitrogen). The cDNA was then amplified using indexing primers. Following quantitative real-time PCR (qRT-PCR) amplification, the library size was determined using the Bioanalyzer (Agilent Technologies 2100) and library concentrations were assessed using qRT-PCR (EvaGreen). The libraries were then processed by Novogene. Single-end reads [50 base-pair (bp) read length] were sequenced using the Illumina Novaseq 6000 sequencer to produce at least 20 M clean reads per sample. The adaptor sequences were trimmed and bases with a Phred quality score less than 20 were removed.

 Prediction of Mature miRNA Sequences

Raw sequencing reads were trimmed to remove adapters using Cutadapt 2.8 software (Martin, 2011). Trimmed reads were processed using the mapper.pl function of miRDeep2 software (Friedländer et al., 2012) to remove reads shorter than 18 nucleotides and collapsed reads. The processed reads were then mapped onto the genome assembly of O. melastigma (Om_v0.7.RACA) with Bowtie 1.2.2 software (Langmead, 2010), using the following parameters: -n 1 -e 80 -l 18 -a –best –strata, producing mapped reads in BWT format. The BWT files were converted to mirDeep2 ARF format and parsed to remove unmatched nucleotides at the 3′ end, using the convert_bowtie_output.pl and parse_mappings.pl functions of mirDeep2 respectively. The parsed mapped reads were used to predict miRNA sequences using the miRDeep2.pl function of mirDeep2. The miRNA prediction results from all samples were merged and parsed into a single list of unique predicted mature miRNA sequences, which were then matched to known miRNAs by nucleotide BLAST to the miRBase database (Kozomara et al., 2018) using the blastn function of BLAST + 2.6.0 (Camacho et al., 2009) with the following parameters: -task blastn-short -strand plus -num_alignments 1. The BLAST hits were filtered using the following parameters: the first 2–7 nucleotides are identical between the predicted miRNA and known miRNA sequences, hit coverage was ≥90%, no indels present, and a maximum of two mismatches.

 Expression Analysis of Small RNAs

Trimmed sequencing reads were remapped onto the mature miRNA sequences using Bowtie with the following parameters: -v 2 -a -m 1 –norc –best –strata, with mapped reads output in SAM format. The SAM files were converted to BAM format using samtools (Li et al., 2009). The number of reads mapped to each mature miRNA sequence was then counted, producing a read counts table. Differential expression analysis of the miRNA read counts was performed using DESeq2 (Love et al., 2014). Differentially expressed miRNAs (DEmiRNAs) were determined using cutoff values of | log2 (fold change: treatment/control)| > 1 and B&H corrected p-value < 0.05.

 miRNA Target Genes Prediction

Target mRNAs of significant differentially expressed miRNAs were predicted using IntaRNA v2:3:1 (Camacho et al., 2009; Li et al., 2009; Langmead, 2010; Martin, 2011; Friedländer et al., 2012; Fagegaltier et al., 2014; Love et al., 2014; Wright et al., 2014; Mann et al., 2017; Busch et al., 2018; Kozomara et al., 2018; Raden et al., 2018)[20, 21] and miRanda v3:3a (John et al., 2005), with default parameters. Target mRNAs were treated as true targets when they met the following criteria: (1) targets were predicted by both tools; (2) interaction sites overlapped between predictions.

 Transcriptome Sequencing and Bioinformatic Analysis

Total RNA extracted from the same set of embryos (three replicates each from normoxic and hypoxic conditions) were subjected to cDNA library construction. The libraries were processed by Novogene. Paired-end reads (150 bp × 2) were sequenced on the Illumina Novaseq 6000 sequencer to produce at least 23 M clean reads per sample. The raw sequencing reads were checked for quality using FastQC (Andrews, 2010). The sequence reads were then mapped onto the genome assembly of O. melastigma (Om_v0.7.RACA) using Spliced Transcripts Alignment to a Reference (STAR) 2.7.1a (Dobin et al., 2013), and aligned reads were output in BAM format. Aligned reads were annotated to exon features of the Om_v0.7.RACA gene annotation file using the htseq-count function of HTseq 0.11.2 (Anders et al., 2015), and the resulting read-count data were subjected to differential expression analysis using DESeq2. Genes with | log2 (fold change: treatment/control)| > 1 and B&H corrected p-value < 0.05 were considered differentially expressed genes (DEGs).

 Assignment of Human Gene Symbols to Differentially Expressed Gene and Functional Analysis of Differentially Expressed Gene

Differentially expressed gene were assigned human gene symbols by reciprocal best hit (RBH) matching of the human and O. melastigma proteomes. Briefly, protein BLAST of the human (GRCh38.p12) protein sequences was performed against the medaka (Om_v0.7.RACA) protein sequences, and vice versa, using BLAST + 2.6.0 (Camacho et al., 2009). The top hit, with the lowest e-value, for each query sequence was extracted and compared between the two sets of BLAST queries; Protein sequences where the query and top hit matched in both directions were identified as RBHs. The DEGs were then subjected to gene ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using DAVID 6.81, enriched biological processes and KEGG pathways (p < 0.05) were identified.

 Complementary DNA Synthesis and Quantitative PCR Analysis

One μg of total RNA was converted to cDNA by using SuperScript™ VILO™ cDNA Synthesis Kit (Invitrogen). Then qPCR analysis was performed on the cDNA from normoxia and hypoxia groups using StepOnePlus™ Real-Time PCR System (Thermofisher). The primer sequencing was shown as Supplementary Table 2.

 Results Exposure to Hypoxia Did Not Affect the Hatchability and Survival Rate of Medaka Embryos

Following normoxia or hypoxia exposure, the viability and hatchability of the embryos were measured for 4 weeks (28 days). Our results indicated that exposure to hypoxia did not affect the hatchability (Figure 1A) and survival rates (Figure 1B) of the medaka embryos, as compared to the normoxic group.

Hypoxia exposure has no effect on the hatchability and viability of medaka embryos. (A) Embryo hatchability was assessed continuously for 28 days after normoxia or hypoxia exposure. (B) Total embryo viability was determined at day 28 after normoxia or hypoxia exposure. n.s. represents not statistically significant (p > 0.05).

 Hypoxia Was Predicted to Alter Embryonic Development Through miRNA Dysregulation

To identify changes to the miRNA profile in the hypoxia-exposed embryos, small RNA sequencing was employed. For this, we obtained 72.8 and 79.3 million quality-trimmed raw reads from normoxic and hypoxic-exposed embryos, respectively, translating to a total of 4.36 Gb of qualified data (Table 1). The clean sequencing reads were mapped to the O. melastigma genome to determine the miRNA content of marine medaka embryos. We identified 253 conserved miRNAs and 572 novel miRNAs (Supplementary Table 3) in the medaka embryo. When we compared the expression level of these miRNAs in the normoxic and hypoxic-exposed embryos, we found 131 miRNAs with altered regulation, including 50 upregulated and 81 downregulated miRNAs (Figure 2A and Table 2). Of these, 75 and 56 dysregulated miRNAs were conserved and novel miRNAs, respectively (Table 2).

Statistics for small RNA sequencing.

Sample name Normoxia 1 Normoxia 2 Normoxia 3 Hypoxia 1 Hypoxia 2 Hypoxia 3
Number of raw reads 26,132,796 25,721,292 21,373,916 22,951,657 34,874,183 21,931,340
Number of trimmed reads 26,059,732 25,589,068 21,114,713 22,877,502 34,752,221 21,650,001
Number of raw bases 1,306,639,800 1,286,064,600 1,068,695,800 1,147,582,850 1,743,709,150 1,096,567,000
Number of trimmed bases 746,846,191 755,768,344 609,584,250 659,360,186 984,978,164 604,696,719

Hypoxia induces miRNA dysregulation, leading to altered medaka embryonic development. (A) Volcano plot indicating differential miRNA expression in medaka embryos after hypoxia or normoxia exposure. miRNA with | log2 (fold change: hypoxia/normoxia)| > 1 and –log B&H corrected p-value > 1.3 were considered differentially expressed. Red dots represent upregulated miRNA, green dots represent downregulated miRNA, and black dots represent miRNA with no significant change. (B) Rich factor plot showing the biological processes altered by the miRNAs target predicted by Gene Ontology (GO) enrichment analysis of Database for Annotation, Visualization and Integrated Discovery (DAVID). Dot size represents the number of target genes. The color intensity of each dot represents the significance of the biological process.

Hypoxia dysregulated miRNA in medaka embryos.

Conserved/novel miRNA miRNA sequence log2 fold change (hypoxia/normoxia) Adjusted p-value
Novel miRNA_409 TGTGTAAAGAGATGGTCGACGGT 4.96 0.0035
Novel miRNA_428 TACATTGGATCTGTGTGTGCGCGCC 4.93 0.016
Novel miRNA_569 TCATGACTGACTACCTGGCTCAGGT 4.93 0.011
Novel miRNA_18 ACCCTCCGGATCGGCTGCCTCCGCC 4.72 0.0065
Novel miRNA_66 TTCGGTTTCCTGTGTTTTACATC 4.29 0.0028
Novel miRNA_503 TGCTCTGAGGGACGTGGACTTGG 4.24 0.0054
Novel miRNA_94 TGCTCCCCTCTCCCTGCTCCTGG 4.17 0.0060
Novel miRNA_518 GATCATGCGGACACAGCGGGCCT 3.97 0.043
Novel miRNA_296 TTCGGTCTGGGGCTGACGCTCA 3.89 0.00018
Novel miRNA_422 TCCGTCGGTCCTGCTCGGGTCC 3.80 4.39E-05
Novel miRNA_349 CTGTCCGTCATTCTGCACTCGC 3.78 1.76E-11
Novel miRNA_565 TGCCTGTGGATCCCTTGTGAAGAG 3.49 1.25E-06
Novel miRNA_258 TGGATTGTGAAGGAGAAAAGCG 3.42 0.043
Novel miRNA_311 AGGACCCGACTTGTAACTTTGA 3.31 1.15E-07
Novel miRNA_380 TTAGCGGTCCACTAACTTTGTTT 3.27 0.028
Novel miRNA_416 TGTGTGTCTGTAGATCAGGAGC 3.18 0.0033
Novel miRNA_338 TCTCCAGATCTGGTCTCTAGTC 2.91 1.58E-05
Novel miRNA_345 TCTGAACTGCAGCGCCACCTGC 2.89 0.040
Novel miRNA_106 TGAGTGGTGCTGTAGCTGGCTGGCC 2.77 0.00012
Novel miRNA_537 TGAAGAGTGTGCAGCGGCTGGACGT 2.77 0.0030
xla-miR-210-5p AAGCCACTGACTAACGCACATT 2.68 1.64E-07
Novel miRNA_319 TTGATGTCTACTTGGGTTCTGGAGT 2.61 0.00067
Novel miRNA_494 GGTCTCTAGTCTCCAGGTCTGG 2.61 0.0050
Novel miRNA_197 TGAGGTTGGGAGCTCAGACGGG 2.53 0.0054
Novel miRNA_158 TCGTGTATCAGATCTAGGGAGACTA 2.33 0.00038
Novel miRNA_137 TTTGAGAAGACTCGGAGGCGGTGG 2.33 0.043
Novel miRNA_423 TTCAGATCTCCTGTAGAGCAG 2.15 0.041
Novel miRNA_145 TCCGTGTCTCTCCTGTGCCTGGCG 2.03 0.0024
Novel miRNA_510 TTTGGTGTGGCCTGCTGTGTGTC 2.02 0.0060
Novel miRNA_507 TTTGTGACCTGTTATACTGCTA 2.02 0.0076
Novel miRNA_320 TTGGATAAACTGCTGTAACTCAGC 2.02 0.0022
Novel miRNA_450 TCAAGGACGTGCTGGTCAAGG 1.96 0.011
Novel miRNA_134 TCGATTGTGGAGGAGAAAAGTG 1.95 0.0090
Novel miRNA_548 TTCTGTAACGAAGAGGCTGAGACGT 1.90 0.0019
Novel miRNA_290 TCAGACTTGTACAGACTCCGCAGGG 1.87 0.00028
Novel miRNA_207 TCCAATCTCAGGGTGAACACTCTA 1.87 0.0065
Novel miRNA_417 TGAAGACTTTTGGGATTATCTAAA 1.86 0.017
Novel miRNA_414 TTGGCGTGATGCGAGCTTGGCTTG 1.83 4.39E-05
Novel miRNA_556 CGAGACCTGGAGTTTAACATCT 1.69 0.035
Novel miRNA_60 TTTCTTGGGTCTTGCATGGACA 1.64 0.00012
Novel miRNA_292 TCCAATCTCAGGGTGGACACTCTA 1.59 0.021
cli-miR-1388-3p ATCTCAGGTTCGTCAGCCCATG 1.49 0.0085
Novel miRNA_62 TAGGATAATGAGGTCTCTTTAAGG 1.43 0.017
gmo-let-7d-5p TGAGGTAGTTGGTTGTATGGTT 1.41 0.0079
Novel miRNA_360 TATTGTTGATTGGTGGAATCCA 1.37 0.016
Novel miRNA_44 ACCCCAACATGTAGCACTTACT 1.37 0.0084
Novel miRNA_210 TGTGACTGAAGCGCTTTGGGCCCTC 1.36 0.0040
Novel miRNA_31 CGTGTAATGTGCGTTCAAGAAC 1.34 0.040
oga-miR-25 CATTGCACTTGTCTCGGTCTGA 1.16 0.017
Novel miRNA_344 AAAGTGCTTCTTGTTGGGTTGG 1.15 0.0079
ocu-miR-183-5p TATGGCACTGGTAGAATTCACT –1.03 0.0079
oga-miR-26b TTCAAGTAATCCAGGATAGGTT –1.18 0.042
Novel miRNA_247 CCGGGAGTGGGACTGTTTGCACT –1.18 0.013
mmr-miR-204 TTCCCTTTGTCATCCTATGCCT –1.18 0.014
gmo-miR-216b-5p TAATCTCTGCAGGCAACTGTGA –1.19 0.042
oga-miR-17 CAAAGTGCTTACAGTGCAGGTAG –1.25 0.041
xla-miR-27b-3p TTCACAGTGGCTAAGTTCTGCA –1.25 0.016
ocu-miR-7a-5p TGGAAGACTAGTGATTTTGTTGTT –1.35 0.00078
xla-miR-203-3p GTGAAATGTTTAGGACCACTTG –1.36 0.0024
oni-miR-217 TACTGCATCAGGAACTGATTGGC –1.36 0.017
oga-miR-16 TAGCAGCACGTAAATATTGGC –1.37 0.0012
xla-miR-200b-3p TAATACTGCCTGGTAATGATGAT –1.38 0.00029
Novel miRNA_139 TTTGAGTGCGGTACCACATCTG –1.41 0.016
Novel miRNA_547 GCTGCTCAGCACTCCAAACGCG –1.50 0.011
Novel miRNA_363 GATTTCAGTGGATTGAAGAGTA –1.54 0.0020
ocu-miR-205-5p TCCTTCATTCCACCGGAGTCTG –1.54 0.0024
ccr-miR-10c TACCCTGTAGATCCGGATTTGTG –1.60 0.0013
gmo-miR-15b-5p TAGCAGCGCATCATGGTTTGAAAC –1.67 0.00027
nle-miR-30c TGTAAACATCCTACACTCTCAGCT –1.69 8.55E-06
cli-miR-455-5p TATGTGCCCTTGGACTACATCGT –1.71 0.00010
hhi-miR-10d TACCCTGTAGAACCGAATGTGTG –1.74 0.00078
mle-miR-1-3p TGGAATGTAAAGAAGTATGTAT –1.78 0.0060
ocu-miR-10b-5p TACCCTGTAGAACCGAATTTGTG –1.80 0.0019
nle-miR-454 TAGTGCAATATTGCTTATAGGGTGT –1.82 0.00044
dre-miR-181a-5-3p ACCATCGACCGTTGACTGTGCC –1.83 0.00016
oga-miR-30d TGTAAACATCCCCGACTGGAAGCT –1.84 1.40E-05
xla-miR-130c-5p GCCCTTTTTCTGTTGCACTACT –1.85 0.0022
gmo-miR-126-3p TCGTACCGTGAGTAATAATGCA –1.94 0.00066
pny-miR-725 TTCAGTCATTGTTTCTGGTCGT –1.97 2.56E-05
ocu-miR-107-3p AGCAGCATTGTACAGGGCTATCA –2.01 0.0084
oga-miR-190 TGATATGTTTGATATATTAGGTTG –2.06 0.00021
ssa-miR-148a-3p TCAGTGCATTACAGAACTTTGTT –2.12 0.0061
gmo-miR-152-3p TCAGTGCATAACAGAACTTTG –2.13 1.39E-05
gmo-miR-20b-5p AAAAGTGCTCACAGTGCAGATA –2.14 0.00029
Novel miRNA_33 TGGGCTCACATCAACCTCTTCATC –2.17 0.00076
nbr-miR-29c ACTGATTTCCTCTGGTGCTTAGA –2.18 0.0077
Novel miRNA_20 GTTTTTTTAGGTTTTGATTTTC –2.23 2.94E-05
pny-miR-7132a-3p TGAGGCGTTTAGAACAAGTTCA –2.26 0.00037
nle-miR-20 TAAAGTGCTTATAGTGCAGGTAG –2.27 1.54E-05
oga-miR-181a AACATTCAACGCTGTCGGTGAGT –2.38 3.68E-06
Novel miRNA_135 CAAACCATAATGTGCTGCCTCT –2.40 0.0069
ipu-miR-219b GGAGTTGTGGATGGACATCACGC –2.41 0.016
cli-miR-181a-2-3p ACCATCGACCGTTGACTGTACC –2.44 2.17E-07
ola-miR-301b-5p GCTCTGACAATGTTGCACTACT –2.46 1.78E-05
Novel miRNA_168 CAGAACTTAGTTCATTAGTGAGCA –2.57 0.0024
ocu-miR-124-5p CGTGTTCACAGCGGACCTTGATT –2.62 8.60E-05
ola-miR-92b TATTGCACTTGTCCCGGCCTCC –2.68 0.0060
ola-miR-194-3p CCAGTGGAGGTGCTGTTACCTG –2.70 0.00010
pny-miR-218b TTGTGCTTGATCTAACCATGCA –2.74 0.033
sbo-miR-214 TACAGCAGGCACAGACAGGCAG –2.77 0.0071
ocu-miR-181b-5p AACATTCATTGCTGTCGGTGGGTT –2.77 1.99E-08
gmo-miR-216a-3p CACAATGGCCTCTGGGATTATG –2.78 4.39E-05
ocu-miR-499-5p TTAAGACTTGCAGTGATGTTTA –2.79 5.14E-05
cpo-miR-200a-3p TAACACTGTCTGGTAACGATGTT –2.80 8.84E-10
novel miRNA_530 AATGAGACTGATGTTCTTCCTCTGC –2.83 1.41E-08
gmo-miR-19d-3p TGTGCAAACCCATGCAAAACTG –2.83 8.95E-06
oni-miR-101b GTACAGTACTATGATAACTGAA –2.83 2.78E-06
pny-miR-458 ATAGCTCTTTAAATGGTACTGC –2.85 1.97E-06
oga-miR-19b TGTGCAAATCCATGCAAAACTG –2.93 4.09E-16
ocu-miR-153-5p TCATTTTTGTGATGTTGCAGCT –2.95 0.011
oga-miR-30e TGTAAACATCCTTGACTGGAAGCT –2.96 1.88E-09
ocu-miR-140-5p CAGTGGTTTTACCCTATGGTAG –3.07 1.41E-08
oga-miR-18 TAAGGTGCATCTAGTGCAGATAG –3.15 2.57E-08
pny-miR-106 TAAAGTGCTTACAGTGCAGGTAG –3.16 1.43E-10
pny-miR-301a TAGTGCAATAGTATTGTCAAAGC –3.18 0.033
gmo-miR-7552-5p TTACAATTAAAGGATATTTCTG –3.19 3.65E-05
gmo-miR-449a-5p AGGCAGTGTCTCGTTAGCTGGA –3.28 0.0066
gmo-miR-93-5p AAAAGTGCTGTTTGTGCAGGTAG –3.30 1.48E-14
gmo-miR-135d-5p TATGGCTTTTTATTCCTACGTGA –3.34 2.46E-05
pny-miR-734 TAAATGCTGCAGAATTGTGCTC –3.50 1.43E-16
ocu-miR-199a-5p CCCAGTGTTCAGACTACCTGTTC –3.52 9.46E-07
ssa-miR-2188-3p GCTGTGTGAGGTCAGACCTATC –3.58 0.00011
ocu-miR-133a-5p AGCTGGTAAAATGGAACCAAATC –3.84 0.027
gmo-miR-181d-5p AACATTCATTGCTGTCGCTGGGTT –3.92 5.83E-06
xla-miR-449-5p AGGCAGTGCAATGTTAGCTGGC –3.96 0.033
abu-miR-135a-5p TATGGCTTTTTATTCCTATCTGA –3.99 2.39E-12
pny-miR-8159 TCAGTAACTGGAATCTGTCCCTGCA –4.06 0.0016
pny-miR-219a AGAATTGTGTATGGACATCTGT –4.43 4.86E-23
Novel miRNA_208 TTACTGTTTTATCTCTTATTTTTAG –4.51 0.00049
gga-miR-205c-5p ACTTCACACCACTGAAATCTGG –5.04 0.013
mle-miR-219-5p TGATTGTCCAAACGCAATTCTTG –5.60 7.09E-13

Bioinformatic analysis using the miRanda and IntaRNA algorithms was used to predict the dysregulated miRNA target genes. By combining the results of the 2 algorithms, we found that the hypoxia-dysregulated miRNAs were predicted to target 2466 genes in the embryo (Supplementary Table 4). The predicted target genes were subjected to Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 analysis to further understand the effect of hypoxia-dysregulated miRNA. Our results showed that the hypoxia exposure could alter the gene cluster involved in many biological processes related to embryo development, such as cartilage development, microtubule cytoskeleton organization, liver development, neural crest cell migration, myelination in the peripheral nervous system, and pectoral fin development (Figure 2B). More importantly, fertility and sex determination-related cell signaling including the retinoic acid receptor signaling pathway and the steroid hormone-mediated signaling pathway were highlighted in our analysis (Figure 2B). Although, the hypoxia exposure had no effect on the GSI and HSI indexes (Supplementary Figure 1).

 Hypoxia Altered Embryonic Development and Sex Determination Through the Regulation of miRNA–mRNA Pairs

As miRNAs are one of the major mediators of gene expression, we applied comparative transcriptomic analysis to determine the differential gene expression occurring as a result of hypoxia exposure. A total of 163.6 million quality-trimmed raw reads, translating to 8.18 Gb of data were obtained from the RNA sequencing (Table 3). In the comparative transcriptome analysis, we found 3009 DEGs, including 1543 upregulated genes and 1466 downregulated genes in the embryos exposed to hypoxia (Figure 3A and Supplementary Table 5).

Statistics for transcriptome sequencing.

Sample name Normoxia 1 Normoxia 2 Normoxia 3 Hypoxia 1 Hypoxia 2 Hypoxia 3
Number of raw reads 25,665,909 30,833,254 31,167,240 25,663,817 26,464,103 23,853,887
Number of bases 7,699,772,700 9,249,976,200 9,350,172,000 7,699,145,100 7,939,230,900 7,156,166,100
% GC content 50.08 49.84 49.66 50.48 50.32 50.61
Number of reads mapped to genome 23,444,759 28,474,120 28,582,198 23,973,604 24,686,360 21,986,391

Hypoxia-induced miRNA–mRNA dysregulation, resulting in developmental defects and alteration of sex determination and differentiation in medaka embryos. (A) Volcano plot showing differential mRNA expression in embryos after hypoxia or normoxia exposure. mRNAs with | log2 (fold change: hypoxia/normoxia)| > 1 and –log B&H corrected p-value > 1.3 were considered differentially expressed. Red dots represent upregulated mRNA, green dots represent downregulated mRNA, and black dots represent mRNA with no significant change. (B) Rich factor plots showing the biological processes altered by the miRNA–mRNA pairs using Gene Ontology (GO) enrichment analysis of Database for Annotation, Visualization and Integrated Discovery (DAVID). Dot size represents the number of genes. The color intensity of each dot represents the significance of the biological process. (C) Hypoxia-dysregulated miRNA–mRNA pairs related to gonad development, sex hormone response, male sex characterization. The complementary sequence shows the interaction between miRNA and mRNA. (D) Quantitative PCR validated the differential gene expression in embryo after the hypoxia exposure. * Represented a statistically significant different of gene expression between normoxia and hypoxia groups.

We then integrated the miRNA and mRNA sequencing results to look at the reversed expression of miRNA and mRNA. Our results indicated that hypoxia could lead to the dysregulation of 167 miRNA–mRNA interaction pairs including 31 downregulated miRNA-upregulated mRNA pairs (Table 4) and 136 upregulated miRNA-downregulated mRNA pairs (Table 5). The miRNA-targeted mRNA was subjected to DAVID analysis to understand the effect of hypoxia-dysregulated miRNA–mRNA pairs in the embryo. Our results showed that hypoxia exposure could result in alterations to miRNA-mediated genes related to different developmental processes, such as adipose tissue development, skeletal muscle tissue growth, and post-embryonic development (Figure 3B). Additionally, hypoxia exposure altered the gene cluster related to neuron migration and amacrine cell differentiation, leading to the interference of brain development through the control of miRNA–mRNA interactions (Figure 3B). Finally, we further assessed the possible impact of hypoxia on sex determination and differentiation through miRNA-mRNA regulation. Using gene ontology analysis, we identified miRNA-mRNA pairs that may be involved in biological functions related to sex determination and differentiation including novel miR-145-col9a3 and novel miRNA-94- arid5b in gonad development, novel miRNA-210-ca2, novel miRNA-106-nr2f2, nbr-miR-29c-nr4a1, and ola-miR-92b-akr1d1 in sex hormone response, and novel miRNA-145-mns1, nle-miR-20-sord, and ipu-miR-219b-abcc8 in sex characteristic development (Figure 3C). The differential gene expression was further validated by using quantitative PCR (qPCR), our data showed that the result of qPCR matched with the finding of RNA sequencing (Figure 3D). Taken together, our data suggest that hypoxia may alter sex determination and differentiation through the regulation of miRNA-mRNA pairs.

Hypoxia-induced downregulated miRNA-upregulated mRNA pairs in medaka embryos.

Downregulated miRNA Upregulated mRNA
ccr-miR-10c rimbp2
gmo-miR-93-5p MPZL3
ocu-miR-199a-5p ap3d1
sbo-miR-214 ngef
sbo-miR-214 atp10b
sbo-miR-214 hdac4
sbo-miR-214 best1
sbo-miR-214 GPR137B
gmo-miR-181d-5p CDKL5
gmo-miR-181d-5p grm1a
oga-miR-26b ahcyl1
ocu-miR-7a-5p si:dkey-288a3.2
abu-miR-135a-5p f2rl1.2
Novel miRNA_247 ablim2
Novel miRNA_247 poln
Novel miRNA_247 bmper
nbr-miR-29c nr4a1
oga-miR-19b CTSS
nle-miR-20 sord
mmr-miR-204 fnip1
mmr-miR-204 hkdc1
mmr-miR-204 tfeb
ipu-miR-219b abcc8
gmo-miR-449a-5p FBLN2
gmo-miR-449a-5p vwc2
gmo-miR-449a-5p kcnd2
ocu-miR-124-5p shtn1
ocu-miR-124-5p efhc2
Novel miRNA_547 ppp1r9alb
ola-miR-92b ahdc1
ola-miR-92b AKR1D1

Hypoxia-induced upregulated miRNA-downregulated mRNA pairs in medaka embryos.

Upregulated miRNA Downregulated mRNA
Novel miRNA_503 mdga1
Novel miRNA_503 dhx57
Novel miRNA_503 ulk1a
Novel miRNA_503 nkain2
Novel miRNA_503 adck1
Novel miRNA_503 chp2
Novel miRNA_349 unm_sa821
Novel miRNA_349 apobec2a
Novel miRNA_349 stk36
Novel miRNA_349 lrig1
Novel miRNA_349 PDZRN4
Novel miRNA_349 neurod1
Novel miRNA_349 cdh15
Novel miRNA_416 bnc2
Novel miRNA_416 mcm3
Novel miRNA_537 VAT1
Novel miRNA_537 CHRNA1
Novel miRNA_537 zgc:63863
Novel miRNA_414 JMY
Novel miRNA_414 tepsin
Novel miRNA_414 c1qbp
Novel miRNA_414 txlnba
Novel miRNA_565 atp6v0a2a
Novel miRNA_565 gmpr2
Novel miRNA_292 dhx36
Novel miRNA_292 WDR77
Novel miRNA_510 acp6
Novel miRNA_510 prmt9
Novel miRNA_510 chodl
Novel miRNA_510 ANXA2
Novel miRNA_510 asmt
Novel miRNA_296 lrp5
Novel miRNA_338 her8.2
Novel miRNA_338 si:dkey-40m6.8
Novel miRNA_338 tp53inp1
Novel miRNA_338 EPHA6
Novel miRNA_94 ccdc15
Novel miRNA_94 myom2a
Novel miRNA_94 sphkap
Novel miRNA_94 ogfrl1
Novel miRNA_94 mep1b
Novel miRNA_94 zmp:0000000760
Novel miRNA_94 lap3
Novel miRNA_94 pprc1
Novel miRNA_94 ttll6
Novel miRNA_94 sema3c
Novel miRNA_94 pcsk2
Novel miRNA_94 map3k12
Novel miRNA_94 s1pr3a
Novel miRNA_94 itga10
Novel miRNA_94 adgra2
Novel miRNA_94 lhx4
Novel miRNA_94 arid5b
Novel miRNA_94 casp2
Novel miRNA_94 clybl
Novel miRNA_94 fkbpl
Novel miRNA_94 gtf2f1
Novel miRNA_94 c1qtnf6b
Novel miRNA_94 tmem204
Novel miRNA_94 rab33a
Novel miRNA_94 uncx4.1
Novel miRNA_94 phf2
Novel miRNA_319 cass4
Novel miRNA_319 MARK1
Novel miRNA_319 adamts14
Novel miRNA_319 SLC4A11
Novel miRNA_319 ppp1r9a
Novel miRNA_290 dis3l2
Novel miRNA_507 klhl2
Novel miRNA_409 actr5
Novel miRNA_409 rbm24a
Novel miRNA_18 ino80
Novel miRNA_18 FRMPD4
Novel miRNA_18 fam234b
Novel miRNA_18 tonsl
Novel miRNA_18 smarcc1b
Novel miRNA_18 ccne2
Novel miRNA_450 crb1
Novel miRNA_134 hic1
Novel miRNA_134 tox2
Novel miRNA_134 lrp3
Novel miRNA_134 cul2
Novel miRNA_134 nup93
Novel miRNA_197 vgll2b
Novel miRNA_197 fam124b
Novel miRNA_197 ino80da
Novel miRNA_197 RANGAP1
Novel miRNA_197 tprkb
Novel miRNA_197 carmil2
Novel miRNA_197 pik3cd
gmo-let-7d-5p ints6l
gmo-let-7d-5p si:dkey-119f1.1
gmo-let-7d-5p lpl
gmo-let-7d-5p nwd2
gmo-let-7d-5p iqsec1b
Novel miRNA_345 ints6
Novel miRNA_345 pik3ap1
Novel miRNA_210 zfr2
Novel miRNA_210 ca2
Novel miRNA_210 glra2
Novel miRNA_210 PLXNA2
Novel miRNA_210 si:dkey-175g6.2
Novel miRNA_518 btbd9
oga-miR-25 msh3
Novel miRNA_344 sox6
Novel miRNA_494 usp28
Novel miRNA_494 slx4
Novel miRNA_494 PIK3CA
Novel miRNA_494 plxnb1a
Novel miRNA_494 aldh18a1
Novel miRNA_106 slc35e4
Novel miRNA_106 pex5lb
Novel miRNA_106 tlx2
Novel miRNA_106 babam1
Novel miRNA_106 PTPRD
Novel miRNA_106 nr2f2
Novel miRNA_106 pcloa
Novel miRNA_106 elavl4
Novel miRNA_106 cacng6b
Novel miRNA_106 ccdc106a
Novel miRNA_422 pcsk1
Novel miRNA_422 tppp2
Novel miRNA_422 zgpat
Novel miRNA_422 gtf2h1
Novel miRNA_145 col9a3
Novel miRNA_145 cdk5rap1
Novel miRNA_145 mns1
Novel miRNA_145 rgmb
Novel miRNA_145 syt4
Novel miRNA_145 gap43
Novel miRNA_145 actn2b
Novel miRNA_145 si:ch211-283g2.1
Novel miRNA_145 neurod4
Novel miRNA_290 si:ch211-87j1.4
Novel miRNA_290 dync1li2
Novel miRNA_290 chrnd
Discussion

Using small-RNA sequencing analysis, we identified 253 conserved miRNA and 572 novel miRNAs in the medaka embryo. The larger number of novel miRNAs being identified in this study compared to conserved miRNAs suggests that many marine medaka miRNAs are yet to be identified. When we compared our results to the limited miRNAs previously identified in marine medaka which is approximately 200 novel miRNAs identified in different organs such as brain, liver, and gonads (Lai et al., 2016). Many novel miRNAs in embryos are important in fish embryonic development, such as bone and gonadal development (Hausser and Zavolan, 2014; Li et al., 2016).

We then looked at the miRNAs altered by hypoxia exposure by comparing the miRNA profile of normoxic and hypoxic exposed embryos. We observed some conserved miRNAs, which are reported to play important roles in embryonic development. For example, the downregulated miR-214 is a developmental regulator, which controls the polycomb protein, Ezh2, in skeletal muscle and embryonic stem cells (Gan et al., 2016). Additionally, an in vitro study demonstrated that miR-214 expression levels are controlled by the transcription factor Twist-1 in the development of specific neural cell populations (Presslauer et al., 2017). Another hypoxia-downregulated miRNA − miR-29c reportedly affects lateral development and cardiac circulation through the Wnt4/β-catenin signaling pathway in zebrafish (Juan et al., 2009). The involvement of miR-29c has also been implicated in bovine blastocyst development and embryo implantation of rats with endometriosis (Lee et al., 2009; Shen et al., 2020), suggesting the importance of hypoxia dysregulated miR-29c in embryonic development.

Our results also highlighted miR-19b, which was reduced following hypoxia exposure, and is reported to impair cardiac development in zebrafish by targeting ctnnb1(Goossens et al., 2013). Our results in medaka are also concordant with a bird study of miRNAs in great tits that found miRNA-19b to be a hypoxia-responsive miRNA, which regulated MAPK1 expression in embryonic fibroblasts (Cai et al., 2018). In addition, the expression level of miRNA-19b is also associated with embryo quality (Li et al., 2014). We also found hypoxia to suppress the expression of miR-204; this miRNA is reported to alter neuronal migration and cortical morphogenesis during embryonic development in mouse embryos (Chen et al., 2018), and both human and zebrafish studies have demonstrated miR-204-mediated control of developmental lymphangiogenesis (Abu-Halima et al., 2020). The reduction of these miRNAs by hypoxia is suggestive of the possible effects of hypoxia exposure on embryonic development through miRNA regulation. Following analysis of miRNAs for which expression was induced by hypoxia exposure, we found elevated expression levels in a large number of novel miRNAs, but not conserved miRNAs. This result made determining the possible effect of the novel miRNA profile change complex. Therefore, we applied two algorithms (MiRanda and IntaRNA) to predict miRNA target genes. To strengthen our findings, we also conducted comparative transcriptome sequencing followed by the miRNA and mRNA integrative analysis, to further determine miRNA-mRNA interaction pairs. Hundred and ninety-nine miRNA–mRNA pairs were identified, and DAVID analysis was performed on the miRNA target genes to determine the effect of hypoxia-dysregulated miRNA-mRNA pairs. For the data analysis, we focused primarily on the endpoint of embryonic development and sex determination and differentiation. In the functional characterization, novel miRNA-145-col9a3 was reported to be involved in gonad development, as Col9a3 encodes the α3 chain of type IX collagen, which is expressed during mouse testis development (Venø et al., 2017; Jung et al., 2019). Another mouse study demonstrated the specified expression of Col9a3 in testicular cords during the early stages of gonadal differentiation (Hanson-Kahn et al., 2018). More importantly, Col9a3 expression was markedly upregulated in the male, but remained very low in the female during embryo development (Hanson-Kahn et al., 2018), suggesting the possible role of Col9a3 in sex differentiation.

Our findings also highlighted some hypoxia-dysregulated miRNA-mRNA pairs such as novel miRNA-210-ca2, novel miRNA-106-nr2f2, and nbr-miR-29c-nr4a1, which could contribute to sex-steroid hormone response. The novel miRNA, miRNA-210, mediates carbonic anhydrase II (CA2), which is a metalloenzyme responsible for the maintenance of the acid-base balance in body systems (McClive and Sinclair, 2003). Rat studies have demonstrated an association between CA2 and sex-steroid hormones, for example, CA2 expressed in the rat lateral prostate and seminal vesicles is under testosterone regulation (Perera et al., 2001), and it is involved in bicarbonate production − a function that particularly characterizes the rat lateral prostate (Perera et al., 2001). The expression level of CA2 is also differentially regulated by testosterone in the dorsal and lateral prostate in rats (Sanyanga et al., 2019). Other than male sex hormones, CA2 is also regulated by female sex hormones. Hepatic carbonic anhydrase is reportedly induced by estrogen in rat models (Härkönen and Väänänen, 1988), and estrogen and progesterone differentially regulate CA2 in ovariectomized rat uteri (Härkönen et al., 1991). The other hypoxia-altered miRNA-mRNA pair, novel miRNA-106-nr2f2, was also found to be associated with the development of sex organs. Nuclear Receptor Subfamily 2 Group F Member 2 (nr2f2), encodes a ligand-inducible transcription factor and is expressed during early ovarian development of embryogenesis in ovarian somatic cells (Garg, 1975). A review of genetic disease studies showed that rare disease differences or disorders of sex development (DSD) are associated with mutations in NR2F2 (Karim et al., 2016). Furthermore, a study of human genetic disorders demonstrated that NR2F2 loss of function caused defects in testis development in individuals with DSD (Rastetter et al., 2014).

As well as novel miRNAs, our result also highlighted the conserved miR-29c-nr4a1, which is also involved in sex differentiation. Nr4a1, an orphan nuclear receptor, is important for hormone-induced steroidogenesis in Leydig cells, where it plays a pivotal role in regulating the expression of several genes involved in male sex differentiation (Kremen and Chan, 2019). Nr4a1 can be activated by many transcription factors such as myocyte enhancer factor 2 to control Leydig cell gene expression (Bashamboo et al., 2018). It has also been reported that Nr4a1 is a regulator of Insulin-like 3 (INSL3) transcription (Garg, 1975), which is a hormone produced by fetal Leydig cells of the testis to regulate testicular descent during fetal life. In adults, Nr4a1 acts as a germ cell survival factor (Abdou et al., 2016).

Finally, novel miR-145-mns1 and nle-miR-20-sord pairs are reported to play roles in male sex characteristics. Meiosis-specific nuclear structural protein 1 (Mns1) is necessary for spermiogenesis and motile cilia function, such as sperm flagella, through its interaction with ciliary proteins (Robert et al., 2006; Daems et al., 2014). Studies of Mns1-deficient mice demonstrate that homozygous loss-of-function mutations in Mns1 resulted in laterality defects and male infertility (Boschen et al., 2018). This finding was supported by a human population study combining genome-wide SNP mapping and whole-exome sequencing analysis that identified MNS1 variant as a cause of male infertility in humans (Zhou et al., 2012).

For nle-miR-20-sord, sorbitol dehydrogenase (Sord) expression is reported to be regulated by androgens in the human prostate, suggesting that this is the location of the physiological role of sord (Ta-Shma et al., 2018). Furthermore, Sord expression is induced to support epididymal sperm motility during spermatogenic differentiation in mouse spermatogenesis (Leslie et al., 2020). Accordingly, reduced expression of Sord in the mouse epididymis results in inhibition of sperm maturation through the impairment of secretory functions of the epididymis (Szabó et al., 2010).

 Conclusion

Our integrative analysis of small RNA sequencing and transcriptome sequencing identified a cluster of miRNA–mRNA pairs that may involve in embryonic development. Our findings have also provided important insights into the molecular mechanisms through which miRNAs interact with different target genes that may regulate the reproductive axis and mediate hypoxia-altered plasticity of sex determination and differentiation in fish. But further studies are required to confirm the role of miRNA–mRNA pairs in physiological processes of hypoxia-altered sex biased in fish.

 Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA706118.

 Ethics Statement

The animal study was reviewed and approved by all animal research procedures were approved by the Animal Ethics Committee on Research Experiments involving Animal Subjects (A-0244) of the City University of Hong Kong.

 Author Contributions

KL, SC, and RK contributed to conception and design of the study. NT, YK, and WT organized the database. CT, CL, and YK performed the experiments. YC, XL, and TC performed the bioinformatics and statistical analysis. KL and RK wrote the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Funding

This study was fully supported by a grant from the General Research Fund (CityU 11102918) of the Research Grants Council of Hong Kong SAR, People’s Republic of China. KL was supported by the Hong Kong SAR, Macao SAR, and Taiwan Province Talented Young Scientist Program of Guangxi.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.736362/full#supplementary-material

 References Abdou H. S. Robert N. M. Tremblay J. J. (2016). Calcium-dependent Nr4a1 expression in mouse Leydig cells requires distinct AP1/CRE and MEF2 elements. J. mol. Endocrinol. 56 151–161. 10.1530/JME-15-0202 26647388 Abu-Halima M. Khaizaran Z. A. Ayesh B. M. Fischer U. Khaizaran S. A. Al-Battah F. (2020). MicroRNAs in combined spent culture media sterility. Fertil. Steril. 113 970–980. 10.1016/j.fertnstert.2019.12.028 32222254 Anders S. Pyl P. T. Huber W. (2015). HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics (Oxford, England). 31 166–169. 10.1093/bioinformatics/btu638 25260700 Andrews S. (2010). A Quality Control Tool for High Throughput Sequence Data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ Baroiller J. F. D’Cotta H. Saillant E. (2009). Environmental effects on fish sex determination and differentiation. Sex Dev. 3 118–135. 10.1159/000223077 19684457 Bashamboo A. Eozenou C. Jorgensen A. Bignon-Topalovic J. Siffroi J. P. Hyon C. (2018). Loss of function of the nuclear receptor NR2F2, encoding COUP-TF2, causes testis development and cardiac defects in 46,XX children. Am. J. Hum. Genet. 102 487–493. 10.1016/j.ajhg.2018.01.021 29478779 Boschen K. E. Gong H. Murdaugh L. B. Parnell S. E. (2018). Knockdown of Mns1 increases susceptibility to craniofacial defects following gastrulation-stage alcohol exposure in mice. Alcohol. Clin. Exp. Res. 42 2136–2143. 10.1111/acer.13876 30129265 Bouchareb A. Le Cam A. Montfort J. Gay S. Nguyen T. Bobe J. (2017). Genome-wide identification of novel ovarian-predominant miRNAs: new insights from the medaka (Oryzias latipes). Sci. Rep. 7:40241. 10.1038/srep40241 28071684 Busch A. Andreas S. A. Backofen R. (2018). IntaRNA: ecient prediction of bacterial sRNA targets incorporating target site accessibility and seed regions. Bioinformatics 24 2849–2856. 10.1093/bioinformatics/btn544 18940824 Cai H. Zhu X. X. Li Z. F. Zhu Y. P. Lang J. H. (2018). MicroRNA dysregulation and steroid hormone receptor expression in uterine tissues of rats with endometriosis during the implantation window. Chin. Med. J. 131 2193–2204. 10.4103/0366-6999.240808 30203794 Camacho C. Coulouris G. Avagyan V. Ma N. Papadopoulos J. Bealer K. (2009). BLAST+: architecture and applications. BMC Bioinformatics 10:421. 10.1186/1471-2105-10-421 20003500 Chen X. Qu Y. Cheng Y. Wang J. Lei X. Song G. (2018). MiR-19b-3p Regulates MAPK1 expression in embryonic fibroblasts from the great tit (parus major) under hypoxic conditions. Cell. Physiol. Biochem. 46 546–560. 10.1159/000488621 29617681 Cheung C. H. Chiu J. M. Wu R. S. (2014). Hypoxia turns genotypic female medaka fish into phenotypic males. Ecotoxicology 23 1260–1269. 10.1007/s10646-014-1269-8 25011919 Daems C. Martin L. J. Brousseau C. Tremblay J. J. (2014). MEF2 is restricted to the male gonad and regulates expression of the orphan nuclear receptor NR4A1. Mol. Endocrinol. (Baltimore, Md.) 28 886–898. 10.1210/me.2013-1407 24694307 Dobin A. Davis C. A. Schlesinger F. Drenkow J. Zaleski C. Jha S. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics (Oxford, England) 29 15–21. 10.1093/bioinformatics/bts635 23104886 Eshel O. Shirak A. Dor L. Band M. Zak T. Markovich-Gordon M. (2014). Identification of male-specific amh duplication, sexually differentially expressed genes and microRNAs at early embryonic development of Nile tilapia (Oreochromis niloticus). BMC Genomics 15:774. 10.1186/1471-2164-15-774 25199625 Fagegaltier D. König A. Gordon A. Lai E. C. Gingeras T. R. Hannon G. J. (2014). Genome-wide survey of sexually dimorphic expression of Drosophila miRNAs identifies the steroid hormone-induced miRNA let-7 as a regulator of sexual identity. Genetics 198 647–668. 10.1534/genetics.114.169268 25081570 Feng Y. P. Chen J. F. Huang P. Wang X. Wang J. Peng X. L. (2014). Expression analysis of differentially expressed miRNAs in male and female chicken embryos. Genet. Mol. Res. 13, 3060–3068. Friedländer M. R. Mackowiak S. D. Li N. Chen W. Rajewsky N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40 37–52. 10.1093/nar/gkr688 21911355 Gan S. Huang Z. Liu N. Su R. Xie G. Zhong B. (2016). MicroRNA-140-5p impairs zebrafish embryonic bone development via targeting BMP-2. FEBS Lett. 590 1438–1446. 10.1002/1873-3468.12190 27130837 Garg L. C. (1975). Induction of hepatic carbonic anhydrase by estrogen. J. Pharmacol. Exp. Ther. 192 297–302. 804034 Gay S. Bugeon J. Bouchareb A. Henry L. Delahaye C. Legeai F. (2018). MiR-202 controls female fecundity by regulating medaka oogenesis. PLoS Genet. 14:e1007593. 10.1371/journal.pgen.1007593 30199527 Goossens K. Mestdagh P. Lefever S. Van Poucke M. Van Zeveren A. Van Soom A. (2013). Regulatory microRNA network identification in bovine blastocyst development. Stem Cells Dev. 1 1907–1920. 10.1089/scd.2012.0708 23398486 Hanson-Kahn A. Li B. Cohn D. H. Nickerson D. A. Bamshad M. J. University of Washington Center for Mendelian Genomics (2018). Autosomal recessive Stickler syndrome resulting from a COL9A3 mutation. Am. J. Med. Genet. 176 2887–2891. 10.1002/ajmg.a.40647 30450842 Härkönen P. L. Mäkelä S. I. Valve E. M. Karhukorpi E. K. Väänänen H. K. (1991). Differential regulation of carbonic anhydrase II by androgen and estrogen in dorsal and lateral prostate of the rat. Endocrinology 128 3219–3232. 10.1210/endo-128-6-3219 1903702 Härkönen P. L. Väänänen H. K. (1988). Androgen regulation of carbonic anhydrase II, a major soluble protein in rat lateral prostate tissue. Biol. Reprod. 38 377–384. 10.1095/biolreprod38.2.377 3129039 Hausser J. Zavolan M. (2014). Identification and consequences of miRNA-target interactions beyond repression of gene expression. Nat. Rev. Genet. 15, 599–612. Jing J. Wu J. J. Liu W. Xiong S. T. Ma W. G. Zhang J. (2014). Sexbiased miRNAs in gonad and their potential roles for testis development in yellow catfish. PLoS One 9:e107946. 10.1371/journal.pone.0107946 25229553 John B. Enright A. J. Aravin A. Tuschl T. Sander C. Marks D. S. (2005). Human MicroRNA targets. PLoS Biol. 2:e363. 10.1371/journal.pbio.0020363 15502875 Juan A. H. Kumar R. M. Marx J. G. Young R. A. Sartorelli V. (2009). Mir-214-dependent regulation of the polycomb protein Ezh2 in skeletal muscle and embryonic stem cells. Mol. Cell 36 61–74. 10.1016/j.molcel.2009.08.008 19818710 Juanchich A. Le Cam A. Montfort J. Guiguen Y. Bobe J. (2013). Identification of differentially expressed miRNAs and their potential targets during fish ovarian development. Biol. Reprod. 88:128. 10.1095/biolreprod.112.105361 23595902 Jung H. M. Hu C. T. Fister A. M. Davis A. E. Castranova D. Pham V. N. (2019). MicroRNA-mediated control of developmental lymphangiogenesis. eLife 8:e46007. 10.7554/eLife.46007 31478836 Kang L. Cui X. Zhang Y. Yang C. Jiang Y. (2004). Identification of miRNAs associated with sexual maturity in chicken ovary by Illumina small RNA deep sequencing. BMC Genomics 14:352. 10.1186/1471-2164-14-352 23705682 Karim K. Giribabu N. Muniandy S. Salleh N. (2016). Estrogen and progesterone differentially regulate carbonic anhydrase II, III, IX, XII, and XIII in ovariectomized rat uteri. Syst. Biol. Reprod. Med. 62 57–68. 10.3109/19396368.2015.1112699 26709452 Kozomara A. Birgaoanu M. Griffiths-Jones S. (2018). Mirbase: from microrna sequences to function. Nucleic Acids Res. 47 D155–D162. 10.1093/nar/gky1141 30423142 Kremen J. Chan Y. M. (2019). Genetic evaluation of disorders of sex development: current practice and novel gene discovery. Curr. Opin. Endocrinol. Diabetes obes. 26 54–59. 10.1097/MED.0000000000000452 30507699 Langmead B. (2010). Aligning short sequencing reads with Bowtie. Curr. Protoc. Bioinformat. 2010:Chater11,Unit–11.7. 10.1002/0471250953.bi1107s32 21154709 Lai K. P. Li J. W. Tse A. C. K. Chan T. F. Wu R. S. S. (2016). Hypoxia alters steroidogenesis in female marine medaka through miRNA regulation. Aquat. Toxicol. 172 1–8. Lau K. Lai K. P. Bao J. Y. Zhang N. Tse A. Tong A. (2014). Identification and expression profiling of microRNAs in the brain, liver and gonads of marine medaka (Oryzias melastigma) and in response to hypoxia. PLos One 9:e110698. 10.1371/journal.pone.0110698 25350659 Lee Y. B. Bantounas I. Lee D. Y. Phylactou L. Caldwell M. A. Uney J. B. (2009). Twist-1 regulates the miR-199a/214 cluster during development. Nucleic Acids Res. 37 123–128. 10.1093/nar/gkn920 19029138 Leslie J. S. Rawlins L. E. Chioza B. A. Olubodun O. R. Salter C. G. Fasham J. (2020). MNS1 variant associated with situs inversus and male infertility. Eur. J. Hum. Genet. 28 50–55. 10.1038/s41431-019-0489-z 31534215 Li H. Handsaker B. Wysoker A. Fennell T. Ruan J. Homer N. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25 2078–2079. 10.1093/bioinformatics/btp352 19505943 Li J. W. Lin X. Tse A. Cheung A. Chan T. F. Kong R. Y. (2016). Discovery and functional characterization of novel miRNAs in the marine medaka Oryzias melastigma. Aquat. Toxicol. (Amsterdam, Netherlands) 175 106–116. 10.1016/j.aquatox.2016.03.013 27002527 Li M. Hu X. Zhu J. Zhu C. Zhu S. Liu X. (2014). Overexpression of miR-19b impairs cardiac development in zebrafish by targeting ctnnb1. Cell. Physiol. Biochem. 33 1988–2002. 10.1159/000362975 25034767 Love M. I. Huber W. Anders S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. 10.1186/s13059-014-0550-8 25516281 Ma H. Hostuttler M. Wei H. R. Rexroad C. E. Yao J. B. (2012). Characterization of the rainbow trout egg microRNA transcriptome. PLoS One 7:e39649. 10.1371/journal.pone.0039649 22761856 Mann M. Wright P. R. Backofen R. (2017). IntaRNA 2.0: enhanced and customizable prediction of RNA-RNA interactions. Nucleic Acids Res. 45 W435–W439. 10.1093/nar/gkx279 28472523 Martin M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.J. 17 10–12. 10.1089/cmb.2017.0096 28715235 McClive P. J. Sinclair A. H. (2003). Type II and type IX collagen transcript isoforms are expressed during mouse testis development. Biol. Reprod. 68 1742–1747. 10.1095/biolreprod.102.008235 12606408 Mishima T. Takizawa T. Luo S. S. Ishibashi O. Kawahigashi Y. Mizuguchi Y. (2008). MicroRNA (miRNA) cloning analysis reveals sex differences in miRNA expression profiles between adult mouse testis and ovary. Reproduction 136, 811–822. Navarro-Martín L. Viñas J. Ribas L. Díaz N. Gutiérrez A. Di Croce L. (2011). DNA methylation of the gonadal aromatase (cyp19a) promoter is involved in temperature-dependent sex ratio shifts in the European sea bass. PLoS Genet. 7:e1002447. 10.1371/journal.pgen.1002447 22242011 Pelley J. (2006). Dead zones might masculinize fish. Environ. Sci. Technol. 40:2862. 10.1021/es0630137 16719079 Perera E. M. Martin H. Seeherunvong T. Kos L. Hughes I. A. Hawkins J. R. (2001). Tescalcin, a novel gene encoding a putative EF-hand Ca(2+)-binding protein, Col9a3, and renin are expressed in the mouse testis during the early stages of gonadal differentiation. Endocrinology 142 455–463. 10.1210/endo.142.1.7882 11145610 Pörtner H. O. Peck M. A. (2010). Climate change effects on fishes and fisheries: towards a cause-and-effect understanding. J. Fish. Biol. 77 1745–1779. 10.1111/j.1095-8649.2010.02783.x 21078088 Presslauer C. Tilahun Bizuayehu T. Kopp M. Fernandes J. M. Babiak I. (2017). Dynamics of miRNA transcriptome during gonadal development of zebrafish. Sci. Rep. 7:43850. 10.1038/srep43850 28262836 Raden M. Ali S. M. Alkhnbashi O. S. Busch A. Costa F. Davis J. A. (2018). Freiburg RNA tools: a central online resource for RNA-focused research and teaching. Nucleic Acids Res. 46 W25–W29. 10.1093/nar/gky329 29788132 Rastetter R. H. Bernard P. Palmer J. S. Chassot A. A. Chen H. Western P. S. (2014). Marker genes identify three somatic cell types in the fetal mouse ovary. Dev. Biol. 394 242–252. 10.1016/j.ydbio.2014.08.013 25158167 Reddon A. R. Hurd P. L. (2013). Water pH during early development influences sex ratio and male morphin a West African cichlid fish, Pelvicachromis pulcher. Zoology (Jena) 116 139–143. 10.1016/j.zool.2012.11.001 23474178 Robert N. M. Martin L. J. Tremblay J. J. (2006). The orphan nuclear receptor NR4A1 regulates insulin-like 3 gene transcription in Leydig cells. Biol. Reprod. 74 322–330. 10.1095/biolreprod.105.044560 16237153 Sanyanga T. A. Nizami B. Bishop ÖT. (2019). Mechanism of action of non-synonymous single nucleotide variations associated with α-carbonic anhydrase ii deficiency. Molecules (Basel, Switzerland). 24 3987. 10.3390/molecules24213987 31690045 Shang E. H. H. Yu R. M. K. Wu R. S. S. (2006). Hypoxia affects sex differentiation and development, leading to a male-dominated population in zebrafish (Danio rerio). Environ. Sci. Technol. 40 3118–3122. 10.1021/es0522579 16719120 Shen Y. Lu H. Chen R. Zhu L. Song G. (2020). MicroRNA-29c aing Wnt4. Mol. Med. Rep. 22 4675–4684. 10.3892/mmr.2020.11584 33173954 Szabó Z. Hämäläinen J. Loikkanen I. Moilanen A. M. Hirvikoski P. Väisänen T. (2010). Sorbitol dehydrogenase expression is regulated by androgens in the human prostate. Oncol. Rep. 23 1233–1239. 10.3892/or_00000755 20372835 Ta-Shma A. Hjeij R. Perles Z. Dougherty G. W. Abu Zahira I. Letteboer S. J. F. (2018). Homozygous loss-of-function mutations in MNS1 cause laterality defects and likely male infertility. PLoS Genet. 14:e1007602. 10.1371/journal.pgen.1007602 30148830 Thomas P. Rahman M. S. (2012). Extensive reproductive disruption, ovarian masculinization and aromatase suppression in Atlantic croaker in the northern Gulf of Mexico hypoxic zone. Proc. Biol. Sci. 279 28–38. 10.1098/rspb.2011.0529 21613294 Thrash J. C. Seitz K. W. Baker B. J. Temperton B. Gillies L. E. Rabalais N. N. (2017). Metabolic roles of uncultivated bacterioplankton lineages in the northern Gulf of Mexico “Dead Zone. mBio 8 e1017–e1017. 10.1128/mBio.01017-17 28900024 Venø M. T. Venø S. T. Rehberg K. van Asperen J. V. Clausen B. H. Holm I. E. (2017). Cortical morphogenesis during embryonic development is regulated by miR-34c and miR-204. Front. Mol. Neurosci. 10:31. 10.3389/fnmol.2017.00031 28232790 Wang W. Liu W. Liu Q. Li B. An L. Hao R. (2016). Coordinated microRNA and messenger RNA expression profiles for understanding sexual dimorphism of gonads and the potential roles of microRNA in the steroidogenesis pathway in Nile tilapia (Oreochromis niloticus). Theriogenology 85 970–978. 10.1016/j.theriogenology.2015.11.006 26719037 Wright P. R. Georg J. Mann M. Sorescu D. A. Richter A. S. Lott S. (2014). CopraRNA and IntaRNA: predicting small RNA targets, networks and interaction domains. Nucleic Acids Res. 42 W119–W123. 10.1093/nar/gku359 24838564 Yamamoto Y. Zhang Y. Sarida M. Hattori R. S. Strüssmann C. A. (2014). Coexistence of genotypic and temperature-dependent sex determination in pejerrey Odontesthes bonariensis. PLoS One 9:e102574. 10.1371/journal.pone.0102574 25036903 Zhou J. Yang F. Leu N. A. Wang P. J. (2012). MNS1 is essential for spermiogenesis and motile ciliary functions in mice. PLoS Genet. 8:e1002516. 10.1371/journal.pgen.1002516 22396656

https://david.ncifcrf.gov