Front. Mar. Sci.Frontiers in Marine ScienceFront. Mar. Sci.2296-7745Frontiers Media S.A.10.3389/fmars.2022.930156Marine ScienceOriginal ResearchSea Urchins in Acute High Temperature and Low Oxygen Environments: The Regulatory Role of microRNAs in Response to Environmental StressHanLingshu1WuYanglei2HaoPengfei2DingBeichen2LiYuanxin2WangWenpei2ZhangXianglei2GaoChuang2WangHeng2WangLuo2ZhangWeijie2ChangYaqing2DingDewen1DingJun2*1School of Marine Sciences, Ningbo University, Ningbo, China2Key Laboratory of Mariculture and Stock Enhancement in North China’s Sea, Ministry of Agriculture and Rural Affairs, Dalian Ocean University, Dalian, China
Edited by: Mark Botton, Fordham University, United States
Reviewed by: Lusheng Xin, Chinese Academy of Fishery Sciences (CAFS), China; Silvia Gomez-Jimenez, Consejo Nacional de Ciencia y Tecnología (CONACYT), Mexico
*Correspondence: Jun Ding, dingjun19731119@hotmail.com
This article was submitted to Marine Biology, a section of the journal Frontiers in Marine Science
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.
Strongylocentrotus intermedius is an economically valuable sea urchin species in China. However, its growth and survival are severely constrained by ocean warming and the hypoxia that often accompanies high water temperatures. MicroRNAs (miRNAs) are important regulators of gene expression in response to environmental change. In this study, high-throughput RNA sequencing was used to investigate changes in miRNA expression in S. intermedius under heat (25°C), hypoxia (2 mg/L O2), and combined heat and hypoxia stresses. Twelve small RNAs libraries were constructed and 17, 14, and 23 differentially expressed miRNAs (DEMs) were identified in the heat, hypoxia, and combined stress groups (P<0.05), respectively. Gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway functional analyses of putative target genes of the DEMs suggested that these miRNAs were important in basal metabolism, apoptosis, oxidative stress, and immune-related pathways. By co-analysis with published transcriptome data, key DEMs (miR-193, miR-184, miR-133, miR-125, miR-2008) and their key target genes (EGF3, ABCB4, CYCL, PAN2, CALN) were identified. Quantitative real-time PCR analysis of the expression of 10 DEMs and their key target genes confirmed the RNA sequencing results. These results provide information on gene expression regulation of the molecular mechanisms underlying the response of S. intermedius to multi-cause environmental stresses.
strongylocentrotus intermediusmicroRNAglobal climate changehigh temperature stresshypoxiaIntroduction
Global warming has profoundly affected the nearshore and marine environments. The resulting environmental changes include increased seawater surface temperature, acidification, increased hypoxic zones, and extreme weather, which seriously threaten the reproduction of marine fishery species and the sustainable development of aquaculture (Blencowe, 2006; Xu et al., 2018; García Molinos, 2020). Among them, water temperature is the most fundamental and widespread ecological (environmental) factor that affects the survival of marine organisms (Pinsky et al., 2019). Increased seawater temperatures have been shown to impact the growth, metabolism, reproduction, and development of marine animals (García-Echauri et al., 2020; Gouda and Agatsuma, 2020; Strøm et al., 2020). Hypoxia caused by high temperatures, nutrient enrichment, water pollution, and high-density farming affects wild and farmed aquatic animals (Vaquer-Sunyer and Duarte, 2008; Breitburg et al., 2018). Dissolved oxygen levels of <2.8 mg/L (equivalent to 2 mL O2/L or 91.4 mM) are considered hypoxic. Very high summer temperatures and depleting aquaculture resources lead to hypoxic conditions that can be fatal or sub-lethal for sea urchins and can cause mass mortality (Riedel et al., 2014). Under hypoxic stress, significant changes in the oxygen consumption rate (Tomi et al., 2011) and the expression of immune and metabolic-related genes were found to occur in sea urchins (Suh et al., 2014; Hao et al., 2022). Therefore, hypoxia is another important factor that limits the survival and growth of sea urchins.
Sea urchins are a model organism for embryological development and a globally crucial marine fishery resource species (Nesbit et al., 2019). Among them, Strongylocentrotus intermedius, which has fast growth and excellent qualities, is China’s main sea urchin culture species, accounting for >90% of the total sea urchin culture in China (Chang et al., 2016). However, S. intermedius is a cold-water species that is sensitive to temperature changes. The high temperatures in the northern summers of 2017 and 2020 led to high mortality rates (up to 80%) of S. intermedius cultured in the Shandong and Liaoning provinces (Ding and Chang, 2020). Furthermore, between 1978 and 2014, the monthly average dissolved oxygen content in the Bohai Sea declined (Shi, 2016). These patterns suggest that environmental factors will seriously threaten the sustainable development of the aquaculture industry. Echinoderms have specific molecular biological mechanisms to respond to stress, such as histone modification, transcription, translation, and post-translational modification (Branco et al., 2013; Huo et al., 2018; Huo et al., 2021).
Various non-coding RNAs are known to be involved in post-transcriptional regulation. Among them, microRNAs (miRNAs) are small non-coding RNAs with about 18-25 nucleotides, which are endogenous regulatory RNAs found mainly in eukaryotes. Most miRNAs are transcribed from DNA sequences into primary miRNA, and processed into precursor RNAs by a series of nucleases to obtain mature miRNAs. The mature miRNAs are assembled into RNA-induced silencing complexes that recognize the target mRNAs through complementary base pairing and instruct the silencing complex to degrade or inhibit the translation of the target mRNAs according to the different degrees of complementation (Fabian et al., 2010; Fu, 2014; Zgheib et al., 2017). MiRNAs play irreplaceable roles in cell proliferation and differentiation, cell apoptosis, gene expression, restoration of homeostasis, and target recognition (Bartel, 2009; Leung and Sharp, 2010; Noman et al., 2017; Tian et al., 2019a; Tian et al., 2019b). Many studies have shown that miRNAs are activated in animals under environmental stress. For example, miR-210-5p and miR-92b-3p were highly activated in sea cucumber in response to the environmental pressures of hypoxia and high temperature (Huo et al., 2021). Sun et al. (2019) reported that miR-122, miR-15b-5p, miR-30b, miR-20a-5p, and miR-7b helped maintain the energy requirements of common carp (Cyprinus carpio) under high temperature stress by regulating glycolysis. In sea cucumbers under salinity stress, miR-2008 and miR-92a were shown to respond by regulating proteins and phospholipids (Tian et al., 2019b). In Chinese shrimp (Fenneropenaeus chinensis), novel-mir-76 and novel-mir-193 responded to high pH stress by regulating lipid metabolism, amino acid metabolism, and carbon metabolism pathways (Li et al., 2019). Multiple environmental factors can affect organisms in natural environments. In a previous study (Hao et al., 2022), we performed transcriptome-wide gene expression profiling by RNA sequencing (RNA-seq) of S. intermedius under short-term high temperature, hypoxia, and combined stresses. We found that exposure to the combined stresses resulted in a two-factor additive effect at the transcriptome level and had a broader effect on immune-related pathways in the S. intermedius than a single stress had (Hao et al., 2022). However, no reports on the involvement of miRNAs in the regulatory mechanism of sea urchins under acute environmental stresses have been published so far.
In this study, we used RNA-seq technology to screen and identify differentially expressed miRNAs (DEMs) under high temperature, low oxygen, and combined stresses. We analyzed the obtained DEMs jointly with published transcriptome data to further explore the molecular mechanisms of S. intermedius in response to multiple environmental stresses. The results provide information for selecting and breeding resistant sea urchins and provide a theoretical basis for healthy sea urchin culture.
MethodExperimental Animals and Treatment
Healthy 1-years-old S. intermedius (average test diameter: 35.74 ± 1.35 mm; average test height: 16.34 ± 1.27 mm; weight: 19.25 ± 0.37 g) were bought from Dalian Haibao Fishery Co., Ltd. (38°91′25′′N” 121°60′25′′E) in Dalian, Liaoning Province, China, then immediately transported to the Key Laboratory of Mariculture and Stock Enhancement in North China’s Sea, Ministry of Agriculture and Rural Affairs, Dalian Ocean University. All the S. intermedius were temporarily raised in the same environment for 14 days before the experiment. The breeding conditions are as follows: temperature, 15 - 16 °C; salinity, 31.22 - 31.36 ppt; pH, 8.15 - 8.25; and oxygen content, 8 - 9 mg/L. Sea urchins were fasted for 48 h before the start of the experiment to expel the intestinal contents.
According to the changes in temperature and dissolved oxygen in the Yellow Sea and the Bohai Sea in the past ten years (Song et al., 2020) and previous experiments (Hao et al., 2022), this study selected 25°C and 2 mg/L oxygen as the conditions of temperature and hypoxia stress. Healthy S. intermedius were randomly cultured in high temperature (HT group, 25 °C, 8 mg/L O2), low oxygen (LO group, 15 °C, 2 mg/L O2), combined stress of heat and hypoxia (HL group, 25 °C, 2 mg/L O2) and control group (NC group, 15 °C, 8 mg/L O2), using 6 S. intermedius per group. The stress experiment was carried out in the automatic temperature control system and the dissolved oxygen control system (Figure 1). The temperature control system’s temperature deviation was less than 0.5 °C, and the oxygen concentration deviation of the dissolved oxygen control system was less than 0.5 mg/L (Huo et al., 2018). Previous experiments showed that individuals died at 12 hours of combined stress, and transcriptome sequencing revealed significant changes in genes related to sea urchin immunity and metabolism at 12 hours of acute stress (Hao et al., 2022). Therefore, we selected 12 h as the time of stress. After 12 h, 6 sea urchins in each group were dissected on ice to absorb the coelomic fluid and centrifuged at 3000 rpm for 10 min at 4 °C. The coelomocyte of S. intermedius were collected and immediately stored in an -80 °C for further experiments.
Temperature dissolved oxygen control device.
RNA Extraction, Small RNA Library Construction, Sequencing, and Annotation
Three coelomocyte samples were selected from every group of sea urchins. The extraction, quality testing, and purification of total RNA were performed according to Wang’s method (Wang et al., 2022). The experimental procedure was performed following the standard steps provided by Illumina, including library preparation and sequencing experiments. The Small RNA sequencing library was prepared using TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA). The Illumina HiSeq 2500 (LC Sciences, Houston, Texas, USA) platform was used to sequence the constructed library, and the reading length was 1×50 bp on a single end. Data analysis was based on this literature (Huo et al., 2021). The obtained sequences with 18-25 nt were used to perform BLAST analysis in miR Base 22.0 database (http://microrna.sanger.ac.uk) to identify miRNAs and examine their variations.
Identification, Target Gene Prediction, and Functional Analysis of Differentially Expressed miRNAs (DEMs)
To identify DEMs between treatment group (HT, LO, and HL) and control group (NC), the expression levels of miRNAs in the whole library were normalized, and then DEMs were screened using Student’s T-test. P<0.05 was considered significant. A heatmap was constructed using TBtools and clustered by row scale. Target genes of miRNAs were predicated by using Target Finder (http://www.targetfinder.org) and were subjected to the enrichment analysis of functions and pathways by GO and KEGG database (P-value<0.05).
Interaction Analysis of miRNA-mRNA
Association analysis of DEM and differentially expressed genes (DEGs) were based on published transcriptome data (Database number: PRJNA776993). All possible positively and negatively correlated miRNA-mRNA pairs were predicted using ACGT101-CORR1.1. Based on the comprehensive analysis of DEGs and DEMs, we screened the negatively associated miRNA-mRNA pairs and constructed interaction networks using Cytoscape 2.8.3 software (http://www.cytoscape.org/).
Verification of Differential Expression of miRNAs
To validate the results of RNA-seq, DEMs and target DEGs were randomly selected for real-time fluorescence quantitative PCR (qRT-PCR). Total RNA extraction and reverse transcription were performed by referring to Han et al. (Han et al., 2021). 18s rRNA and U6 RNA were used as internal reference gene for qRT-PCR, which was conducted with the LightCycler96 Real-time System (Roch, Switzerland) and followed by the manufacturer’s instructions of the FastStart Essential DNA Green Master (Roch, Switzerland). The PCR primers were designed and synthesized by Shanghai Sangon Biotech and the primer sequences are shown in Table 1. The qRT-PCR was performed in a 20 μL reaction sample containing 2 μL of cDNA, 10 μL of FastStart Essential DNA Green Master, 6 μL of ddH2O water, and 1 μL (10 mM) of each primer. The reaction conditions were followed by 40 cycles of 95°C for 30 s, 95°C for 5 s, 60°C for 32 s, 95°C for 15 s, 60°C for 60 s, 95°C for 15 s and 60°C for 15 s. Compared with the control gene, fold-change of expression levers was determined using the 2−ΔΔCT method (Livak and Schmittgen, 2001), both in miRNA and mRNA PCR amplification. Statistical analysis was performed using SPSS software version 19.0 (IBM, Armonk, NY, USA). The data are expressed as the mean ± standard error of the mean (SEM) (n=3). In the results of gene expressions, significant differences (P<0.05) for each variable were first detected using the one-way ANOVA test between different groups, followed by Tukey’s HSD test.
Primers used for DEMs and target DEGs verification.
ResultsMiRNA Profiling of S. intermedius Exposed to Heat, Hypoxia, and Combined Stresses
We constructed small RNA libraries from S. intermedius and sequenced them on an Illumina HiSeq2500 platform to identify the miRNAs involved in heat, hypoxia, and combined stresses. Twelve libraries were established from the four experimental groups. A total of 12,049,632 ± 1,609,813 (NC), 9,169,509 ± 1,182,352 (HT), 12,077,449 ± 2,083,876 (LO), and 12,214,860 ± 1,271,035 (HL) raw reads were obtained. After removing low-quality reads, we mapped the clean reads against the Rfam database, which contains rRNA, tRNA, snRNA, snoRNA families, and Repbase to filter unwanted sequences and obtain valid small RNA data. A total of 7,846,735 ± 1,822,577 (NC), 5,549,774 ± 3,009,597 (HT), 8,869,277 ± 1,415,005 (LO), and 9,613,111 ± 885,961 (HL) valid reads were obtained (Table 2). The length distribution statistics for the total number of filtered valid reads showed that most of the reads were 22-nt (Supplementary Figure 1).
Overview of reads from raw data to cleaned sequences.
Sample
Raw reads
3ADT&length filter
Junk reads
Rfam
Repeats
valid reads
rRNA
tRNA
snoRNA
snRNA
other Rfam RNA
NT1
Total
10243015
2605706
17509
1427085
69370
6190556
434776
925238
8042
12264
46765
% of Total
100.00
25.44
0.17
13.93
0.68
60.44
4.24
9.03
0.08
0.12
0.46
uniq
1461697
704280
5788
23737
489
727811
10184
10654
350
571
1978
% of uniq
100.00
48.18
0.40
1.62
0.03
49.79
0.70
0.73
0.02
0.04
0.14
NT2
Total
12573998
3257425
23380
1739864
102857
7550284
663282
969640
10924
24813
71205
% of Total
100.00
25.91
0.19
13.84
0.82
60.05
5.28
7.71
0.09
0.20
0.57
uniq
1807890
910655
6569
27639
538
862934
12084
11638
479
916
2522
% of uniq
100.00
50.37
0.36
1.53
0.03
47.73
0.67
0.64
0.03
0.05
0.14
NT3
Total
13331885
2220139
25179
1281716
68369
9799364
467576
695791
11079
30746
76524
% of Total
100.00
16.65
0.19
9.61
0.51
73.50
3.51
5.22
0.08
0.23
0.57
uniq
1644836
741590
7730
27300
757
868005
13359
10293
423
697
2528
% of uniq
100.00
45.09
0.47
1.66
0.05
52.77
0.81
0.63
0.03
0.04
0.15
HT1
Total
8577093
1805181
26938
618997
19315
6122514
208024
367369
7720
6963
28921
% of Total
100.00
21.05
0.31
7.22
0.23
71.38
2.43
4.28
0.09
0.08
0.34
uniq
1305399
588644
8583
17730
416
690351
9160
6538
352
335
1345
% of uniq
100.00
45.09
0.66
1.36
0.03
52.88
0.70
0.50
0.03
0.03
0.10
HT2
Total
10530958
1745541
27742
522899
14469
8231847
208230
276875
6876
5066
25852
% of Total
100.00
16.58
0.26
4.97
0.14
78.17
1.98
2.63
0.07
0.05
0.25
uniq
1326233
580640
8908
17655
372
718952
9379
6361
354
285
1276
% of uniq
100.00
43.78
0.67
1.33
0.03
54.21
0.71
0.48
0.03
0.02
0.10
HT3
Total
8400476
5809150
10165
285059
6147
2294962
69102
201759
2146
2271
9781
% of Total
100.00
69.15
0.12
3.39
0.07
27.32
0.82
2.40
0.03
0.03
0.12
uniq
1289897
904019
4195
8924
216
372717
4346
3644
141
130
663
% of uniq
100.00
70.08
0.33
0.69
0.02
28.90
0.34
0.28
0.01
0.01
0.05
LO1
Total
11247361
1550143
17851
897787
28924
8777912
223260
629471
10909
9322
24825
% of Total
100.00
13.78
0.16
7.98
0.26
78.04
1.98
5.60
0.10
0.08
0.22
uniq
1243041
562511
7384
17347
409
655710
8393
6794
401
363
1396
% of uniq
100.00
45.25
0.59
1.40
0.03
52.75
0.68
0.55
0.03
0.03
0.11
LO2
Total
14448447
3232933
14158
872015
25638
10327751
207538
621207
9651
9950
23669
% of Total
100.00
22.38
0.10
6.04
0.18
71.48
1.44
4.30
0.07
0.07
0.16
uniq
1575774
847933
6686
16326
270
704773
7697
6721
384
327
1197
% of uniq
100.00
53.81
0.42
1.04
0.02
44.73
0.49
0.43
0.02
0.02
0.08
LO3
Total
10536539
2202009
17808
811268
24484
7502169
185791
579121
13056
9991
23309
% of Total
100.00
20.90
0.17
7.70
0.23
71.20
1.76
5.50
0.12
0.09
0.22
uniq
1255042
646936
6667
16185
388
585162
7340
6838
431
330
1246
% of uniq
100.00
51.55
0.53
1.29
0.03
46.62
0.58
0.54
0.03
0.03
0.10
HL1
Total
12959469
2435730
14469
531578
19186
9973246
211701
277533
6153
6048
30143
% of Total
100.00
18.79
0.11
4.10
0.15
76.96
1.63
2.14
0.05
0.05
0.23
uniq
1101567
657103
3343
11648
325
429398
6158
3314
294
241
1641
% of uniq
100.00
59.65
0.30
1.06
0.03
38.98
0.56
0.30
0.03
0.02
0.15
HL2
Total
12937862
2080346
16434
572686
32073
10262292
274021
244267
6286
10411
37701
% of Total
100.00
16.08
0.13
4.43
0.25
79.32
2.12
1.89
0.05
0.08
0.29
uniq
950584
505095
3261
12367
364
429768
6608
3394
310
379
1676
% of uniq
100.00
53.14
0.34
1.30
0.04
45.21
0.70
0.36
0.03
0.04
0.18
HL3
Total
10747248
1613140
11916
513863
22544
8603794
215763
261339
4950
6795
25016
% of Total
100.00
15.01
0.11
4.78
0.21
80.06
2.01
2.43
0.05
0.06
0.23
uniq
818902
424313
3156
10891
341
380460
5934
3038
250
253
1416
% of uniq
100.00
51.81
0.39
1.33
0.04
46.46
0.72
0.37
0.03
0.03
0.17
The bioinformatic analysis identified 140 miRNAs in the four experimental groups; 70 known and 70 novel miRNAs (Table 3). By co-expression analysis of the four groups (three comparison groups), a total of 70 (55 co-expressed, HT 1, NC 14), 78 (61 co-expressed, LO 9, NC 8) and 82 (67 co-expressed, HL 13, NC 2) miRNAs were identified in the HT vs NC, LO vs NC, and HL vs NC comparisons, respectively (Figure 2).
Number of known and novel miRNA identified in each sample.
NT libaray
NT1
NT2
NT3
Groups:
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
gp1
2
4
2
4
2
4
gp2a
1
1
1
1
2
2
gp2b
5
4
3
2
4
3
gp3
48
56
49
57
49
57
gp4
12
12
20
19
22
21
HT libaray
HT1
HT2
HT3
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
gp1
2
4
2
4
2
3
gp2a
1
1
1
1
1
1
gp2b
3
2
3
2
1
1
gp3
45
52
43
50
41
48
gp4
21
19
22
21
5
4
LO libaray
LO1
LO2
LO3
Groups:
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
gp1
2
4
2
4
2
4
gp2a
1
1
1
1
1
1
gp2b
4
3
1
1
5
4
gp3
45
52
44
51
44
51
gp4
22
20
27
25
33
28
HL libaray
HL1
HL2
HL3
Groups:
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
Pre-miRNA
Unique miRNA
gp1
2
4
2
4
2
4
gp2a
2
2
1
1
2
2
gp2b
4
3
5
4
4
3
gp3
48
57
48
57
48
56
gp4
26
22
30
25
22
20
Venn of detected miRNA in four groups. (A) HT vs NT, (B) LO vs NT, (C) HL vs NT.
Differentially Expressed miRNAs (DEMs) in Three Comparison Groups
We identified 17, 14, and 23 DEMs in the HT vs NC, LO vs NC, and HL vs NC comparisons, respectively, by analyzing the RNA-seq data. Among the 17 DEMs in the HT vs NC comparison, 7 were up-regulated and 10 were down-regulated. Among the 14 DEMs in the LO vs NC comparison, 8 were up-regulated and 6 were down-regulated. Among the 23 DEMS in the HL vs NC comparison, 14 were up-regulated and 9 were down-regulated (Supplementary Table 1). The heat map (Figure 3) shows that 42 of the DEMs in the three comparisons were significantly differentially expressed (P <0.05), and four of them (PC-3p-56283_77, lva-miR-2008-5p, lva-miR-219-5p_R-1, lva-miR-9-5p_R+1) were significantly different in all three comparison groups.
Differentially expressed miRNAs (DEMs) in the four groups. (A) HT vs NT, (B) LO vs NT, (C) HL vs NT, (D) DEMs (non-repetitive) in the pairwise comparison among the three treatments (P < 0.05). All red rectangles indicate higher levels of miRNAs, and blue rectangles indicate lower levels of miRNAs. Color scale bar represents log2FoldChange.
Identification and Functional Annotation of the Target Genes of the DEMs
The roles of the identified DEMs in regulating the expression and function of their target genes were predicted by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway functional enrichment analysis. A total of 69,912 target genes were predicted for the 17 DEMs identified in the HT vs NC comparison, and 12,825 of them were annotated to 9,879 GO terms under the three main GO categories: biological process (BP), cellular component (CC), and molecular function (MF). Under BP, biological process, signal transduction, and positive regulation of transcription by RNA polymerase II were the three most enriched terms (P <0.05). Under CC, integral component of membrane, membrane, and plasma membrane were the three most enriched terms. Under MF, calcium ion binding, serine-type endopeptidase activity, and GTP binding were the three most enriched terms (Figure 4A and Supplementary Table 2). In addition, 5,710 target genes were annotated to 388 KEGG pathways. Five of these pathways were significantly enriched, namely Glycosaminoglycan biosynthesis - chondroitin sulfate/dermatan sulfate (ko00532), Neuroactive ligand–receptor interaction (ko04080), mTOR signaling pathway (ko04150), FoxO signaling pathway (ko04068), and Drug metabolism - cytochrome P450 (ko00982) (Figure 5A and Supplementary Table 3).
GO enrichment of target genes of DEMs. The x-axis is the rich factor, which indicates the proportion of genes in total genes in a GO term. The y-axis is the gene functional classification of GO. Various colors of plots indicate different values of −log 10 (P-value). Plot diameter represents target gene numbers in a GO term. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.). (A) HT VS NT, (B) LO vs NT, (C) HL vs NT.
KEGG enrichment of target genes of DEMs. The x-axis is the rich factor, which means that the proportion of target genes in total genes in a KEGG term. The y-axis is the gene functional classification of KEGG. Various colors of plots indicate different values of −log 10 (P-value). Plot diameter represents target gene numbers in a KEGG term. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article). (A) HT VS NT, (B) LO vs NT, (C) HL vs NT.
A total of 56,715 target genes were predicted for the 14 DEMs identified in the LO vs NC comparison, and 12,217 of them were annotated to 9,592 GO terms. Under BP and CC, the three most enriched terms were the same as those for the HT vs NC comparison. Under MF, metal ion binding, transferase activity, and nucleic acid binding were the three most enriched terms (Figure 4B). A total of 5,471 target genes were annotated to 383 KEGG pathways. Four of these pathways were the most significantly enriched, namely Histidine metabolism, mTOR signaling pathway, ECM-receptor interaction, and FoxO signaling pathway (Figure 5B).
A total of 79,746 target genes were predicted for the 23 DEMs identified in the HL vs NC comparison, and 14,343 of them were annotated to 9,958 GO terms. Under BP, biological process, oxidation-reduction process, and signal transduction were the three most enriched terms (P <0.05). Under CC, integral component of membrane, plasma membrane, and mitochondrion were the three most enriched terms. Under MF, metal ion binding, zinc ion binding, and transferase activity were the three most enriched terms (P <0.05) (Figure 4C). Unlike the KEGG pathway annotations of the single stresses, the KEGG pathway annotations of the target genes of the 23 DEMs obtained for the combined stresses were mostly enriched in metabolic pathways, such as Glycosaminoglycan biosynthesis, Lysine degradation, Peroxisome, and Valine, leucine and isoleucine degradation and Fatty acid elongation (P <0.05) (Figure 5C).
Prediction of Interactions Between DEMs and Differentially Expressed Genes (DEGs)
The pre-constructed transcriptome library and previously identified DEGs (Hao et al., 2022) were used for co-analysis with the target genes of the DEMs. We identified 15 DEMs that interacted with 33 DEGs under high temperature stress by association analysis. Among them, the expression of 24 DEM–DEG pairs was negatively correlated and the expression of 32 DEM–DEG pairs was positively correlated. We also found 14 DEMs that interacted with 171 DEGs under hypoxia stress; 139 were negatively correlated DEM–DEG pairs and 108 were positively correlated DEM–DEG pairs. Under the combined stresses, 23 DEMs interacted with 102 DEGs; 69 were negatively correlated DEM–DEG pairs and 79 were positively correlated DEM–DEG pairs (Supplementary Table 4). The predicted DEM–DEG network map is shown in Figure 6. Meanwhile, we found that the key DEM in high-temperature stress was lva-miR-193-3p_R-2, whose predicted target genes were DHX33 and CSTF1; the key DEM in hypoxic stress was spu-miR-92e_1ss20TG, with nine essential target genes predicted to be negatively regulated, namely ZNFX1, Pola2, DDB_ G0276821, ACSM3, ctps1, PLOD1, GRID2IP, Pole, and FAU; the key DEM was lva-miR-9- in combined stresses, whose had 28 key negatively regulated target genes, mainly SOD, PAN2, and PRPF.
The network of DEMs and their target DEGs was constructed using Cytoscape 3.3.0. Hexagon is miRNA, circle is target DEGs, red means up-regulated expression, blue means down-regulated expression.
Validation of DEMs and Their Target DEGs by qRT-PCR
A qRT-PCR analysis was performed to validate the RNA-seq results. A total of 10 DEGs (including ABCB4, CYCL, EGF3, ZNFX1, CALN) and 10 DEMs (such as miR-193, miR-184, miR-133, miR-125, miR-2008) were selected for the qRT-PCRs. The results of the qRT-PCR analysis were consistent with the results obtained by the RNA-seq analysis and showed similar expression trends (Figure 7). In addition, the relative expressions of ABCB4, KIF4, ZNFX1, EGF3, DHX33, and PTRPF were significantly increased in high-temperature stress, and CALN, CUL4A, and PAN2 were up-regulated in hypoxic stress; CYCL genes showed a decreasing trend in high temperature and hypoxic stress, but their expressions were increased under combined stress.
qRT-PCR verification of DEMs (A) and target DEGs (B). Each vertical bar represents the Mean ± SD (n=3), U6 and18s rRNA were used as a reference miRNA/gene. *Significant differences at P<0.05 vs control (NT). **Highly significant differences at P<0.01 vs control (NT). Letters above the bars indicate significant differences at P<0.05.
Discussion
Many aquatic organisms are located in environmentally sensitive coastal zones or estuaries, which are extremely vulnerable to climate change. Water temperature and dissolved oxygen are the most important environmental factors that affect the survival of S. intermedius. MiRNAs are known to modulate gene expression and have critical roles in many biological processes, including cell proliferation, apoptosis, differentiation, cell cycle progression, and organ development (Tse et al., 2016; Biggar and Storey, 2018). In this study, miRNA libraries of sea urchins under heat, hypoxia, and combined stresses were constructed, sequenced by RNA-seq, and analyzed together with published transcriptome data to investigate the molecular mechanisms of sea urchins in response to different environmental conditions. The miRNAs were mainly distributed at 22-nt long, which is consistent with a previous study of S. intermedius (Zhan et al., 2018). A total of 140 miRNAs were identified; 70 were known and 70 were newly identified miRNAs. We identified 17, 14, and 23 DEMs in the HT, LO, and HL comparisons with NC, respectively, and found that they were involved in multiple key biological processes related to biosynthesis, metabolism, immunity, and signaling transduction by functional enrichment analysis.
These biological processes may be associated with major changes in the predicted target genes. We found changes in lysosomal pathways related to translation that were similar to the results of Huo et al. (2018). Regarding signal transduction, for the DEMs in the HT vs NC comparison, the mTOR and Wnt signaling pathways were among the enriched pathways; for the DEMs in the LO vs NC comparison, the FoxO signaling and Fanconi anemia pathways were among the enriched pathways; and for the DEMs in the HL vs NC comparison, the Hippo signaling and TGF-β signaling pathways were among the enriched pathways. Interestingly, the KEGG pathway annotations Lysine degradation, Peroxisome, Glycerolipid metabolism, and TGF-β signaling pathway were assigned to targent genes of DEMs in all three comparisons. Small amounts of lysine ingested by animals have been shown to cause triglyceride accumulation and significantly inhibit growth and development (Forster and Ogata, 1998; Ahmed and Khan, 2004). Peroxisomes are involved in multiple metabolic processes, including fatty acid oxidation, ether lipid synthesis, and reactive oxygen species (ROS) metabolism. Recent studies suggest that peroxisomes are critical mediators of cellular responses to various forms of stress, including oxidative stress, hypoxia, starvation, cold exposure, and noise (He et al., 2021a). Jain et al. discovered regulators of cell fitness in high and low oxygen conditions which led to the identification of an essential role for peroxisomes in metabolic adaptation to hypoxic stress (Jain et al., 2020). TGF-β signaling is a multifunctional pathway that controls cell proliferation, cell differentiation and tissue homeostasis (Hata and Chen, 2016). A study has found that TGF-β sensu stricto signaling plays an essential role in the formation of the embryonic skeleton of this sea urchin (Sun and Ettensohn, 2017). This finding together with the changes in the multi-metabolic and apoptosis-multi-species pathways, suggest that the metabolic system of sea urchins will be affected and apoptosis will occur under environmental stresses, in addition, multiple signal transduction pathways may work together to activate the same defense response. The DEMs were also involved in biosynthesis and metabolism-related pathways such as “Glucose metabolism”, “Serotonin metabolism”, and “Galactose metabolism”. These results indicates that sea urchins respond to the adverse effects of environmental stress through a variety of miRNA regulatory mechanisms.
MiRNA sequences are usually partially or totally matched with specific regions of the target genes. This interaction results in endonuclease digestion or translation inhibition, thereby negatively regulating the expression of the target gene. Therefore, identifying target genes that are differentially expressed is critical for understanding the biological functions of the encoded proteins.
Effects of Heat Stress on miRNA Regulation in S. intermedius
MiR-193 plays an important role in cell growth, proliferation, and apoptosis. MiR-193 down-regulates the insulin growth factor 2 gene (IGF2) to induce cell proliferation and migration, which affect the angiogenic process (Yi et al., 2017). We found that miR-193 expression was up-regulated and its target genes DHX33 (putative ATP-dependent RNA helicase DHX33) and CSTF1 (cleavage stimulation factor subunit 1) were down-regulated under high temperature stress. DHX33 is a nucleoprotein involved in cell cycle regulation (Zhang et al., 2011), and CSTF1 encodes a nuclear protein that contains a ribonucleoprotein (RNP)-type RNA binding domain at its N-terminal end, which is essential for mRNA cleavage and polyadenylation. When DNA-damaged cells undergo oxidative stress, CSTF1 can bind to the BARD1/BRCA1 ubiquitin ligase heterodimer, thereby inhibiting 3′ end processing (Fontana et al., 2017). This finding and the regulation of alternative polyadenylation in development, differentiation, and neuronal activation suggest that 3′ end processing can be regulated in response to physiological and pathological stimuli. Therefore, we hypothesized that DNA was damaged under high temperature stress and miR-193a negatively regulated DHX33 and CSTF1 to promote apoptosis.
The expression of miRNA let-7 was also found to be down-regulated under high temperature stress. Members of the let-7 family are significantly elevated under high temperature stimulation, and are involved in regulating cell growth, differentiation, apoptosis, and metabolism (Gibadulinova et al., 2020; Tristán-Ramos et al., 2020; Huo et al., 2021). In the current study, the expression of lva-let-7-p3 was significantly down-regulated under heat stress and combined heat and hypoxic stresses. Lva-let-7-p3 may regulate potential target genes, including PI16 (peptidase inhibitor 16), ATP6V1A (V-type proton ATPase catalytic subunit A isoform X4), and Nup107 (nuclear pore complex protein Nup107 isoform X2).
Effects of Hypoxia on miRNA Regulation in S. intermedius
MiR-184 was shown to be an important regulator of stem cell proliferation and growth (Liu et al., 2015), and its overexpression led to apoptosis and its suppression led to an increase in cell numbers (Foley et al., 2010; Liu et al., 2010). Previous studies have shown that miR-184 overexpression inhibited autophagy and exacerbated oxidative damage, and miR-184 also negatively regulated Wnt signaling in vivo and in vitro (Takahashi et al., 2015). MiR-184 was found to inhibit gene expression in human trabecular meshwork cell cytotoxicity, apoptosis, and the extracellular matrix by targeting the hypoxia-inducible factor HIF-1αin vivo, and it also exhibited angiostatic properties by regulating signaling pathways such as the Akt, TNF-α, and VEGF signaling pathways (Park et al., 2017). Furthermore, reduced expression of miR-184 was shown to inhibit cell growth through the CDC25A-dependent Notch signaling pathway (Cao et al., 2020).
MiR-133 is a key regulator of muscle proliferation and myocardial differentiation and is associated with cell proliferation and apoptosis (Chen et al., 2006; Uchida et al., 2013). The target of miR-133 was shown to be epidermal growth factor (EGF), and when EGF was down-regulated by miR-133, it inhibited its downstream signaling pathways, including the MAPK and AKT signaling pathways (Hernes et al., 2004; Dimauro et al., 2014). MiR-133 was also found to be an important component of the apoptotic pathway. In this study, we found that miR-133 expression was down-regulated under hypoxic stress and that its target genes HSP70, EGF, PAN2, and ABCC9 were up-regulated. In addition, the functional enrichment analysis of miRNA target genes under hypoxic stress showed that most of the genes were related to basal metabolism. Therefore, we hypothesized that miR-133 inhibited apoptosis by regulating HSP70, EGF, and other genes when sea urchins were exposed to hypoxia stress, and that metabolism-related genes were mobilized to maintain vital metabolism processes.
We found that the expression of lva-miR-125-5p was remarkably up-regulated in sea urchins under hypoxic conditions. Previous studies of miR-125-5p focused mainly on its role in regulating cell migration, apoptosis, immunity, proliferation, and cancer (Fassan et al., 2013; Natalia et al., 2018; Xu et al., 2019; Zheng et al., 2019; He et al., 2021b). In zebrafish (Danio rerio), the main target gene for hypoxia response was HIF-1a, which up-regulated the expression of miR-125 under hypoxic conditions (He et al., 2017). However, no up-regulation of HIF genes was detected in the present study. Previous studies have demonstrated that under hypoxic stress, the expression of HIF in the brain of the Wuchang bream (Megalobrama amblycephala) did not differ significantly from that of the control group (Shen et al., 2010). It is speculated that different species may have different hypoxia tolerance thresholds or species-specific differences in oxygen requirements. It might be duo to the short duration of the stress treatment in this study, we plan to further investigate the response of sea urchins under long-term hypoxic stress in future work. MiR-125 was shown to inhibit the expression of the mitofusin Mfn1 and reduce the disordered growth of pulmonary arterial smooth muscle cells under hypoxic conditions, as well as protect pulmonary blood vessels from mitochondrial dysfunction and abnormal remodeling (Ma et al., 2017). These findings provided a theoretical foundation for successful lung organ care. In this study, the potential target genes that may be regulated by lva-miR-125-5p included ZFP708 (zinc finger protein 708-like), ALDH (aldehyde dehydrogenase), KMT2A (lysine methyltransferase 2A), and MAP3K7 (mitogen-activated protein kinase kinase kinase 7). In the Gene Expression Omnibus (GEO) Profiles database (Tanya et al., 2007), ZFP708, ALDH, and MAP3K7 are recorded as being involved in the regulation of organisms in a hypoxic environment. These findings suggest that these genes may be involved in the response of sea urchins to low oxygen environments.
Effects of the Combined Stresses on miRNA Regulation in S. intermedius
Under the combined stresses, spu-miR-2002-3p, which was highly expressed under hypoxic stress, was significantly up-regulated. This miRNA has not been reported in this context in previous studies. We found that under hypoxic stress, 14, 30, 5, 25, 39, and 11 target genes of spu-miR-2002-3p were enriched in neuroactive ligand–receptor interaction, lysine degradation, biosynthesis of unsaturated fatty acids, valine, leucine and isoleucine degradation, lysosome, and tryptophan metabolism KEGG pathways, respectively. The neuroactive ligand–receptor interaction signaling pathway involves all the receptors and ligands on the plasma membrane associated with intracellular and extracellular signaling pathways (Lauss et al., 2007). Under hypoxic stress, the cells pass through the surface of the sea urchin, then the receptors interact with extracellular ligands to trigger a series of metabolic changes. These processes may support the response of sea urchins to high temperature and low oxygen stress to protect them from damage.
We also identified miR-2008 and miR-9, which were significantly expressed under all three conditions. MiR-2008 was identified as a regulator of sea cucumber skin ulcer syndrome outbreaks by deep sequencing, and its target gene is the toll like receptor TLR3 (Li et al., 2012; Zhou et al., 2018). The expression of miR-2008 was also found to be up-regulated in sea cucumber under hypoxic stress. In the present study, miR-2008 expression was up-regulated under all three stress conditions, and the key target genes were EGF3, ABCC9, TUBA, and ACAD. We hypothesized that when S. intermedius is damaged by environmental stresses, the expression of genes related to apoptosis and immunity is significant. Our results support the view that miRNA–target gene interactions are complex, and a single miRNA can regulate multiple target genes simultaneously, and a gene can be regulated by multiple miRNAs simultaneously (Agarwal et al., 2015; Lan et al., 2016; Lai et al., 2022; Miao et al., 2022). The mechanism of the S. intermedius response to environmental stress is also complex. Compared with the effects of a single-factor stress, under multi-factor stress conditions, S. intermedius suffers an increased degree of organismal damage and responds comprehensively through metabolic and immune pathways. The qRT-PCR results further imply that there may be meaningful targeted regulatory interactions between the candidate DEMs and DEGs. In future research, the targeted regulatory interactions of the putative DEM–DEG pairs need to be confirmed by additional in vivo and in vitro experiments (e.g., dual luciferase reporter analysis, gain or loss of function analysis, and western blotting).
Conclusions
In this study, we report the miRNAs profiles of the sea urchin S. intermedius under high temperature, low oxygen, and combined stresses, and identified 17, 14, and 23 DEMs, respectively. By co-analysis with published transcriptome data, key DEMs (miR-193, miR-184, miR-133, miR-125, miR-2008) and their key target genes (EGF3, ABCC9, TUBA, PAN2, CALN) were identified. We found that the S. intermedius defense responses were activated by multiple miRNAs that regulate multiple signal transduction pathways in response to the adverse effects of environmental stress. Furthermore, under the combined heat and hypoxic stresses, the effects of both these factors were superimposed. The expression of target genes associated with apoptosis and oxidative stress was up-regulated, implying transcriptional regulation is a comprehensive response in sea urchins in the face of global climate change. Our results provide information on gene expression regulation of the molecular mechanisms underlying the response of S. intermedius to multi-cause environmental stress and provide a theoretical basis for healthy sea urchin reproduction.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The data presented in the study are deposited in the NCBI repository, accession number PRJNA828018.
Author Contributions
JD, DD, and YC: conceptualization and resources. JD and BD: conceived and designed the experiment. BD, PH, and XZ: performed the experiment. PH, YL, and WW: data curation. LH, PH and YW: analyzed the data. WZ, HW, LW and CG: contributed reagents/materials/analysis tools. YW and LH: writing—original draft preparation. LH: writing—review and editing. JD: funding acquisition. All authors have read and agreed to the published version of the manuscript.
Funding
This work was financially supported by Liaoning Department of Natural Resources Promotion of Marine Economic Development Project No.84[2021], High-level Talent Support Grant for Innovation in Dalian [2020RD03], and Liaoning Province Higher Education Innovation Team and Innovative Talent Support Program Project [LT2019003].
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.
Acknowledgments
We thank Margaret Biswas, PhD, from Liwen Bianji (Edanz) (www.liwenbianji.cn/) for editing the English text of a draft of this manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.930156/full#supplementary-material
Abbreviations
MiRNA, MicroRNA; DEGs, Different expression genes; DEMs, Different expression microRNA; GO, Gene Ontology; BP, Biological processes; MF, Molecular functions; CC, Cellular components; KEGG, Kyoto Encyclopedia of Genes and Genomes; qRT-PCR; Quantitative Real-time PCR.
ReferencesAgarwalV.BellG. W.NamJ. W.BartelD. P. (2015). Predicting Effective microRNA Target Sites in Mammalian mRNAs. 4, e05005. doi: 10.7554/eLife.05005AhmedI.KhanM. A. (2004). Dietary Lysine Requirement of Fingerling Indian Major Carp, Cirrhinus Mrigala (Hamilton). 235, 499–511. doi: 10.1016/j.aquaculture.2003.12.009BartelD. P. (2009). MicroRNAs: Target Recognition and Regulatory Functions. 136, 215–233. doi: 10.1016/j.cell.2009.01.002BiggarK. K.StoreyK. B. (2018). Functional Impact of microRNA Regulation in Models of Extreme Stress Adaptation. 10, 93–101. doi: 10.1093/jmcb/mjx053BlencoweB. J. (2006). Alternative Splicing: New Insights From Global Analyses. 126, 37–47. doi: 10.1016/j.cell.2006.06.023BrancoP. C.BorgesJ. C.SantosM. F.Jensch JuniorB. E.Da SilvaJ. R. (2013). The Impact of Rising Sea Temperature on Innate Immune Parameters in the Tropical Subtidal Sea Urchin Lytechinus Variegatus and the Intertidal Sea Urchin Echinometra Lucunter. 92, 95–101. doi: 10.1016/j.marenvres.2013.09.005BreitburgD.LevinL. A.OschliesA.GrégoireM.ChavezF. P.ConleyD. J.. (2018). Declining Oxygen in the Global Ocean and Coastal Waters. 359, eaam7240. doi: 10.1126/science.aam7240CaoQ.XuW.ChenW.PengD.LiuQ.DongJ.. (2020). MicroRNA-184 Negatively Regulates Corneal Epithelial Wound Healing via Targeting CDC25A, CARM1, and LASP1. 7, 35. doi: 10.1186/s40662-020-00202-6ChangY.TianX.ZhangW.HanF.ChenS.ZhouM.. (2016). Family Growth and Survival Response to Two Simulated Water Temperature Environments in the Sea Urchin Strongylocentrotus Intermedius. 17, 1356. doi: 10.3390/ijms17091356ChenJ. F.MandelE. M.ThomsonJ. M.WuQ.CallisT. E.HammondS. M.. (2006). The Role of microRNA-1 and microRNA-133 in Skeletal Muscle Proliferation and Differentiation. 38, 228–233. doi: 10.1038/ng1725DimauroI.GrassoL.FittipaldiS.FantiniC.MercatelliN.RaccaS.. (2014). Platelet-Rich Plasma and Skeletal Muscle Healing: A Molecular Analysis of the Early Phases of the Regeneration Process in an Experimental Animal Model. 9, e102993. doi: 10.1371/journal.pone.0102993DingJ.ChangY. Q. (2020). Research Progress in Conservation and Utilization of Economic Echinoderm: A Review. 35, 645–656. (in Chinese)FabianM. R.SonenbergN.FilipowiczW. (2010). Regulation of mRNA Translation and Stability by microRNAs. 79, 351–379. doi: 10.1146/annurev-biochem-060308-103103FassanM.PizziM.RealdonS.BalistreriM.GuzzardoV.ZagonelV.. (2013). The HER2-Mir125a5p/Mir125b Loop in Gastric and Esophageal Carcinogenesis. 44, 1804–1810. doi: 10.1016/j.humpath.2013.01.023FoleyN. H.BrayI. M.TivnanA.BryanK.MurphyD. M.BuckleyP. G.. (2010). MicroRNA-184 Inhibits Neuroblastoma Cell Survival Through Targeting the Serine/Threonine Kinase AKT2. 9, 83. doi: 10.1186/1476-4598-9-83FontanaG. A.RigamontiA.LenzkenS. C.FilosaG.AlvarezR.CalogeroR.. (2017). Oxidative Stress Controls the Choice of Alternative Last Exons via a Brahma-BRCA1-CstF Pathway. 45, 902–914. doi: 10.1093/nar/gkw780ForsterI.OgataH. Y. (1998). Lysine Requirement of Juvenile Japanese Flounder Paralichthys Olivaceus and Juvenile Red Sea Bream Pagrus Major. 161, 131–142. doi: 10.1016/S0044-8486(97)00263-9FuX. D. (2014). Non-Coding RNA: A New Frontier in Regulatory Biology. 1, 190–204. doi: 10.1093/nsr/nwu008García-EchauriL. L.LigginsG.Cetina-HerediaP.RoughanM.ColemanM. A.JeffsA. (2020). Future Ocean Temperature Impacting the Survival Prospects of Post-Larval Spiny Lobsters. 156, 104918. doi: 10.1016/j.marenvres.2020.104918García MolinosJ. (2020). Global Marine Warming in a New Dimension. 4, 16–17. doi: 10.1038/s41559-019-1037-5GibadulinovaA.BullovaP.StrnadH.PohlodekK.JurkovicovaD.TakacovaM.. (2020). CAIX-Mediated Control of LIN28/let-7 Axis Contributes to Metabolic Adaptation of Breast Cancer Cells to Hypoxia. 21, 4299. doi: 10.3390/ijms21124299GoudaH.AgatsumaY. (2020). Effect of High Temperature on Gametogenesis of the Sea Urchin Strongylocentrotus Intermedius in the Sea of Japan, Northern Hokkaido, Japan. 525, 151324. doi: 10.1016/j.jembe.2020.151324HanL.QuanZ.HanB.DingB.DingJ. (2021). Molecular Characterization and Expression of the SiUCP2 Gene in Sea Urchin Strongylocentrotus Intermedius. 39, 1523–1537. doi: 10.1007/s00343-020-0181-8HaoP.DingB.HanL.XieJ.WuY.JinX.. (2022). Gene Expression Patterns of Sea Urchins (Strongylocentrotus Intermedius) Exposed to Different Combinations of Temperature and Hypoxia. 41, 100953. doi: 10.1016/j.cbd.2021.100953HataA.ChenY. G. (2016). TGF-β Signaling From Receptors to Smads. 8, a022061. doi: 10.1101/cshperspect.a022061HeA. Y.DeanJ. M.LodhiI. J. (2021a). Peroxisomes as Cellular Adaptors to Metabolic and Environmental Stress. 31, 656–670. doi: 10.1016/j.tcb.2021.02.005HeY.HuangC. X.ChenN.WuM.HuangY.LiuH.. (2017). The Zebrafish miR-125c is Induced Under Hypoxic Stress via Hypoxia-Inducible Factor 1α and Functions in Cellular Adaptations and Embryogenesis. 8, 73846–73859. doi: 10.18632/oncotarget.17994HernesE.FossåS. D.BernerA.OtnesB.NeslandJ. M. (2004). Expression of the Epidermal Growth Factor Receptor Family in Prostate Carcinoma Before and During Androgen-Independence. 90, 449–454. doi: 10.1038/sj.bjc.6601536HeW.ZhangN.LinZ. (2021b). MicroRNA-125a-5p Modulates Macrophage Polarization by Targeting E26 Transformation-Specific Variant 6 Gene During Orthodontic Tooth Movement. 124, 105060. doi: 10.1016/j.archoralbio.2021.105060HuoD.SunL.RuX.ZhangL.LinC.LiuS.. (2018). Impact of Hypoxia Stress on the Physiological Responses of Sea Cucumber Apostichopus Japonicus: Respiration, Digestion, Immunity and Oxidative Damage. 6, e4651. doi: 10.7717/peerj.4651HuoD.SunL.SunJ.ZhangL.LiuS.SuF.. (2021). Sea Cucumbers in a High Temperature and Low Dissolved Oxygen World: Roles of miRNAs in the Regulation of Environmental Stresses. 268, 115509. doi: 10.1016/j.envpol.2020.115509JainI. H.CalvoS. E.MarkhardA. L.SkinnerO. S.ToT. L.AstT.. (2020). Genetic Screen for Cell Fitness in High or Low Oxygen Highlights Mitochondrial and Lipid Metabolism. 181, 716–727. doi: 10.1016/j.cell.2020.03.029LaiK. P.TamN. Y. K.ChenY.LeungC. T.LinX.TsangC. F.. (2022). miRNA–mRNA Integrative Analysis Reveals the Roles of miRNAs in Hypoxia-Altered Embryonic Development- and Sex Determination-Related Genes of Medaka Fish. 8. doi: 10.3389/fmars.2021.736362LanC.ChenQ.LiJ. (2016). Grouping miRNAs of Similar Functions via Weighted Information Content of Gene Ontology. 17, 507. doi: 10.1186/s12859-016-1367-0LaussM.KriegnerA.VierlingerK.NoehammerC. (2007). Characterization of the Drugged Human Genome. 8, 1063–1073. doi: 10.2217/14622416.8.8.1063LeungA. K.SharpP. A. (2010). MicroRNA Functions in Stress Responses. 40, 205–215. doi: 10.1016/j.molcel.2010.09.027LiC.FengW.QiuL.XiaC.SuX.JinC.. (2012). Characterization of Skin Ulceration Syndrome Associated microRNAs in Sea Cucumber Apostichopus Japonicus by Deep Sequencing. 33, 436–441. doi: 10.1016/j.fsi.2012.04.013LiuX.FuB.ChenD.HongQ.CuiJ.LiJ.. (2015). miR-184 and miR-150 Promote Renal Glomerular Mesangial Cell Aging by Targeting Rab1a and Rab31. 336, 192–203. doi: 10.1016/j.yexcr.2015.07.006LiuC.TengZ. Q.SantistevanN. J.SzulwachK. E.GuoW.JinP.. (2010). Epigenetic Regulation of miR-184 by MBD1 Governs Neural Stem Cell Proliferation and Differentiation. 6, 433–444. doi: 10.1016/j.stem.2010.02.017LivakK. J.SchmittgenT. D. (2001). Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2(-Delta Delta C(T)) Method. 25, 402–408. doi: 10.1006/meth.2001.1262LiZ.WangJ.HeY.HuS.LiJ. (2019). Comprehensive Identification and Profiling of Chinese Shrimp (Fenneropenaeus Chinensis) microRNAs in Response to High pH Stress Using Hiseq2000 Sequencing. 50, 14269. doi: 10.1111/are.14269MaC.ZhangC.MaM.ZhangL.ZhangL.ZhangF.. (2017). MiR-125a Regulates Mitochondrial Homeostasis Through Targeting Mitofusin 1 to Control Hypoxic Pulmonary Vascular Remodeling. 95, 977–993. doi: 10.1007/s00109-017-1541-5MiaoB.-B.NiuS.-F.WuR.-X.LiangZ.-B.ZhaiY. (2022). The MicroRNAs-Transcription Factors-mRNA Regulatory Network Plays an Important Role in Resistance to Cold Stress in the Pearl Gentian Grouper. 8. doi: 10.3389/fmars.2021.824533NataliaM. A.AlejandroG. T.VirginiaT. J.Alvarez-SalasL. M. (2018). MARK1 is a Novel Target for miR-125a-5p: Implications for Cell Migration in Cervical Tumor Cells. 7, 54–61. doi: 10.2174/2211536606666171024160244NesbitK. T.FlemingT.BatzelG.PouvA.RosenblattH. D.PaceD. A.. (2019). The Painted Sea Urchin, Lytechinus Pictus, as a Genetically-Enabled Developmental Model. 150, 105–123. doi: 10.1016/bs.mcb.2018.11.010NomanA.FahadS.AqeelM.AliU.AmanullahAnwarS.. (2017). miRNAs: Major Modulators for Crop Growth and Development Under Abiotic Stresses. 39, 685–700. doi: 10.1007/s10529-017-2302-9ParkJ. K.PengH.YangW.KatsnelsonJ.VolpertO.LavkerR. M. (2017). miR-184 Exhibits Angiostatic Properties via Regulation of Akt and VEGF Signaling Pathways. 31, 256–265. doi: 10.1096/fj.201600746RPinskyM. L.EikesetA. M.MccauleyD. J.PayneJ. L.SundayJ. M. (2019). Greater Vulnerability to Warming of Marine Versus Terrestrial Ectotherms. 569, 108–111. doi: 10.1038/s41586-019-1132-4RiedelB.PadosT.PretterebnerK.SchiemerL.SteckbauerA.HaselmairA.. (2014). Effect of Hypoxia and Anoxia on Invertebrate Behaviour: Ecological Perspectives From Species to Community Level. 11, 1491–1518. doi: 10.5194/bg-11-1491-2014ShenR. J.JiangX. Y.PuJ. W.ZouS. M. (2010). HIF-1alpha and -2alpha Genes in a Hypoxia-Sensitive Teleost Species Megalobrama Amblycephala: cDNA Cloning, Expression and Different Responses to Hypoxia. 157, 273–280. doi: 10.1016/j.cbpb.2010.06.013ShiQ. (2016). Spatio-Temporal Mode for Inter-Annual Change of Dissolved Oxygen and Apparent Oxygen Utilization in Summer Bohai Sea. 35, 243–255. (in Chinese)SongG.ZhaoL.ChaiF.LiuF.LiM.XieH. X. (2020). Summertime Oxygen Depletion and Acidification in Bohai Sea, China. 7. doi: 10.3389/fmars.2020.00252StrømJ. F.ThorstadE. B.RikardsenA. H. (2020). Thermal Habitat of Adult Atlantic Salmon Salmo Salar in a Warming Ocean. 96, 327–336. doi: 10.1111/jfb.14187SuhS. S.HwangJ.ParkM.ParkS. Y.RyuT. K.LeeS.. (2014). Hypoxia-Modulated Gene Expression Profiling in Sea Urchin (Strongylocentrotus Nudus) Immune Cells. 109, 63–69. doi: 10.1016/j.ecoenv.2014.08.011SunZ.EttensohnC. A. (2017). TGF-β Sensu Stricto Signaling Regulates Skeletal Morphogenesis in the Sea Urchin Embryo. 421, 149–160. doi: 10.1016/j.ydbio.2016.12.007SunJ.LiuQ.ZhaoL.CuiC.WuH.LiaoL.. (2019). Potential Regulation by miRNAs on Glucose Metabolism in Liver of Common Carp (Cyprinus Carpio) at Different Temperatures. 32, 100628. doi: 10.1016/j.cbd.2019.100628TakahashiY.ChenQ.RajalaR. V. S.MaJ. X. (2015). MicroRNA-184 Modulates Canonical Wnt Signaling Through the Regulation of Frizzled-7 Expression in the Retina With Ischemia-Induced Neovascularization. 589, 1143–1149. doi: 10.1016/j.febslet.2015.03.010TanyaB.TroupD.WilhiteS. E.PierreL.DmitryR.CarlosE.. (2007). NCBI GEO: Mining Tens of Millions of Expression Profiles—Database and Tools Update. 35, 760–765. doi: 10.1093/nar/gkl887TianY.LiG.BuX.ShenJ.TaoZ.ChenL.. (2019a). Changes in Morphology and miRNAs Expression in Small Intestines of Shaoxing Ducks in Response to High Temperature. 46, 3843–3856. doi: 10.1007/s11033-019-04827-2TianY.ShangY.GuoR.ChangY.JiangY. (2019b). Salinity Stress-Induced Differentially Expressed miRNAs and Target Genes in Sea Cucumbers Apostichopus Japonicus. 24, 719–733. doi: 10.1007/s12192-019-00996-yTomiS.StankoviS.LucuČ. (2011). Oxygen Consumption Rate and Na+/K+-ATPase Activity in Early Developmental Stages of the Sea Urchin Paracentrotus Lividus Lam. 65, 431–434. doi: 10.1007/s10152-011-0261-4Tristán-RamosP.Rubio-RoldanA.PerisG.SánchezL.Amador-CuberoS.ViolletS.. (2020). The Tumor Suppressor microRNA Let-7 Inhibits Human LINE-1 Retrotransposition. 11, 5712. doi: 10.1038/s41467-020-19430-4TseA. C. K.Rubio-RoldanA.LiJ. W.WangS. Y.ChanT. F.LaiK. P.. (2016). Hypoxia alters testicular functions of marine medaka through microRNAs regulation. 180, 266–273. doi: 10.1016/j.aquatox.2016.10.007UchidaY.ChiyomaruT.EnokidaH.KawakamiK.TataranoS.KawaharaK.. (2013). MiR-133a Induces Apoptosis Through Direct Regulation of GSTP1 in Bladder Cancer Cell Lines. 31, 115–123. doi: 10.1016/j.urolonc.2010.09.017Vaquer-SunyerR.DuarteC. M. (2008). Thresholds of Hypoxia for Marine Biodiversity. 105, 15452–15457. doi: 10.1073/pnas.0803833105WangH.ZhaoW.LiuX.GangD.ZuoR.HanL.. (2022). Comparative Lipidome and Transcriptome Provide Novel Insight Into Polyunsaturated Fatty Acids Metabolism of the Sea Urchin. 9. doi: 10.3389/fmars.2022.777341XuX.LaiY.ZhouW.HuaZ. (2019). Lentiviral Delivery of a shRNA Sequence Analogous to miR-4319/miR-125-5p Induces Apoptosis in NSCLC Cells by Arresting G2/M Phase. 120, 14017–14027. doi: 10.1002/jcb.28676XuY.RamanathanV.VictorD. G. (2018). Global Warming Will Happen Faster Than We Think. 564, 30–32. doi: 10.1038/d41586-018-07586-5YiF.ShangY. G.LiB.DaiS. L.WuW.ChengL.. (2017). MicroRNA-193-5p Modulates Angiogenesis Through IGF2 in Type 2 Diabetic Cardiomyopathy. 491, 876–882. doi: 10.1016/j.bbrc.2017.07.108ZgheibC.HodgesM. M.HuJ.LiechtyK. W.XuJ. (2017). Long non-Coding RNA Lethe Regulates Hyperglycemia-Induced Reactive Oxygen Species Production in Macrophages. 12, e0177453. doi: 10.1371/journal.pone.0177453ZhangY.ForysJ. T.MiceliA. P.GwinnA. S.WeberJ. D. (2011). Identification of DHX33 as a Mediator of rRNA Synthesis and Cell Growth. 31, 4676–4691. doi: 10.1128/MCB.05832-11ZhanY.LiY.CuiD.PeiQ.SunJ.ZhangW.. (2018). Identification and Characterization of microRNAs From the Tube Foot in the Sea Urchin Strongylocentrotus Intermedius. 4, e00668. doi: 10.1016/j.heliyon.2018.e00668ZhengX.WuZ.XuK.QiuY.SuX.ZhangZ.. (2019). Interfering Histone Deacetylase 4 Inhibits the Proliferation of Vascular Smooth Muscle Cells via Regulating MEG3/miR-125a-5p/IRF1. 13, 41–49. doi: 10.1080/19336918.2018.1506653ZhouX.ChangY.ZhanY.WangX.LinK. (2018). Integrative mRNA-miRNA Interaction Analysis Associate With Immune Response of Sea Cucumber Apostichopus Japonicus Based on Transcriptome Database. 72, 69–76. doi: 10.1016/j.fsi.2017.10.031