Expression Regulation Mechanisms of Sea Urchin (Strongylocentrotus intermedius) Under the High Temperature: New Evidence for the miRNA-mRNA Interaction Involvement

In the context of global warming and continuous high temperatures in the northern part of China during summer, the mortality rate of our main breeding species, Strongylocentrotus intermedius, reached 80% in 2020. How sea urchins respond to high temperatures is of great concern to academia and industry. In this study, we examined the antioxidant enzyme activities of different color tube-footed sea urchins under heat stress and compared their transcriptome and microRNA (miRNA) profiles using RNA-Seq. The results showed that the antioxidant enzyme activities of sea urchins were altered by thermal stress, and the changes in peroxidase activities of red tube-footed sea urchins were particularly significant. Investigations revealed that 1,079 differentially expressed genes (DEGs), 11 DE miRNAs, and 104 “DE miRNA-DEG” pairs in total were detected in sea urchins under high temperature stress. Several mRNA and miRNAs were significantly changed (e.g. HSP70, DnaJ11, HYAL, CALR, miR-184-p5, miR-92a, miR-92c, and miR-124-p5), suggesting these genes and miRNAs exerted important functions in response to high temperature. At the transcriptional level, red tube-footed sea urchins were found to be more sensitive to high temperature and could respond to high temperature rapidly. DE miRNA-mRNA network showed that miR-92b-3p and PC-5p-7420 were the most corresponding miRNAs. Five mRNAs (DnaJ11, SAR1B, CALR, HYOU1, TUBA) may be potential markers of sea urchin response to high temperature. Possible interaction between miRNA-mRNA could be linked to protein folding in the endoplasmic reticulum, Phagosomes, and calcium transport. This study provides a theoretical basis for the molecular mechanism of sea urchin heat tolerance and information that will aid in the selection and breeding of sea urchins with high temperature tolerance.

In the context of global warming and continuous high temperatures in the northern part of China during summer, the mortality rate of our main breeding species, Strongylocentrotus intermedius, reached 80% in 2020. How sea urchins respond to high temperatures is of great concern to academia and industry. In this study, we examined the antioxidant enzyme activities of different color tube-footed sea urchins under heat stress and compared their transcriptome and microRNA (miRNA) profiles using RNA-Seq. The results showed that the antioxidant enzyme activities of sea urchins were altered by thermal stress, and the changes in peroxidase activities of red tube-footed sea urchins were particularly significant. Investigations revealed that 1,079 differentially expressed genes (DEGs), 11 DE miRNAs, and 104 "DE miRNA-DEG" pairs in total were detected in sea urchins under high temperature stress. Several mRNA and miRNAs were significantly changed (e.g. HSP70, DnaJ11, HYAL, CALR, miR-184-p5, miR-92a, miR-92c, and miR-124-p5), suggesting these genes and miRNAs exerted important functions in response to high temperature. At the transcriptional level, red tube-footed sea urchins were found to be more sensitive to high temperature and could respond to high temperature rapidly. DE miRNA-mRNA network showed that miR-92b-3p and PC-5p-7420 were the most corresponding miRNAs. Five mRNAs (DnaJ11, SAR1B, CALR, HYOU1, TUBA) may be potential markers of sea urchin response to high temperature. Possible interaction between miRNA-mRNA could be linked to protein folding in the endoplasmic reticulum, Phagosomes, and calcium transport. This study provides a theoretical basis for the molecular mechanism of sea urchin heat tolerance and information that will aid in the selection and breeding of sea urchins with high temperature tolerance.

INTRODUCTION
Global warming is profoundly affecting the nearshore and marine environments, resulting in environmental changes such as increased seawater surface temperature, acidification and hypoxia, which seriously threaten the population reproduction of marine organisms and the sustainable development of aquaculture (Blencowe, 2006;Sanford and Kelly, 2011;Somero, 2012). Temperature is the most fundamental environmental factor that is extensively affecting the survival of marine life (Harvell et al., 2002). Studies have shown that rising temperature can cause significant changes in the metabolism of marine organisms, which in turn cause changes in their growth, reproduction, behavior, and distribution (García-Echauri et al., 2020;Gouda and Agatsuma, 2020;Strøm et al., 2020). Negative impacts on marine organisms are frequently reported when the environmental temperature increases beyond their thermal optimum (Chan et al., 2019). In recent years, abnormally high temperatures along the coast have caused large-scale deaths of marine organisms, which have become the focus of fishery resource protection and the aquaculture industry (Yasuhara and Danovaro, 2016;Brown et al., 2019).
The sea urchin is a representative species of echinoderms and a model organism for embryonic development, and it is also an important marine fishery resource worldwide (Paredes, 2016). Strongylocentrotus intermedius is the major species of sea urchin cultured in China. However, it is a cold-water species and is sensitive to changes in temperature. The optimal water temperature for S. intermedius in its natural habitat is <20°C, but mortality can occur when the water temperature rises above 25°C (Chang et al., 2016). In recent years, water temperatures in China's northern Yellow and Bohai Seas have frequently exceeded 25°C in summer due to global ocean warming (Zeng et al., 2006). A death rate as high as 80% of S. intermedius cultured in the sea areas adjacent to Shandong and Liaoning resulted from high temperatures in northern China during the summers of 2017 and 2020 (Chang, 2004). The tube foot of the sea urchin is the major tissue that comes in contact with the external environment as it participates in functions such as attachment, movement, sensing, and feeding (Siikavuopio et al., 2012). First discovered in 2010 that the pigment cell composition of S. intermedius tube feet may be split into very distinct red and white groups (Zhang et al., 2010). Thus far, temperature studies on red and white tube-footed sea urchin populations have mostly focused on growth and physiological and biochemical indicators (Zhao et al., 2011;Ding, 2014;Ding et al., 2015;Ding et al., 2017), whereas the molecular mechanisms underlying the differences in thermal resistance between red and white tube-footed sea urchins have not been illuminated.
In recent years, the RNA-Seq method has been widely adopted to investigate the gene response of aquatic animals to temperature stress (Samanta et al., 2006;Jiang et al., 2014;Wang et al., 2019). Studies over the past 5 years have provided important information on the gene expressions of echinoderms that were exposed to temperature stress, and the RNA-Seq results have been used to identify a number of candidate genes associated with the immune response, apoptosis pathway, and metabolic pathways, such as those of the heat shock proteins (HSPs), calmodulin-1 (Cam1), cyclophilins (Cyps) (Gonzalez et al., 2017;Li et al., 2019;Shi et al., 2020;Li et al., 2021). MicroRNAs (miRNAs) are thecriticaly regulators of post-transcriptional gene expression, and they perform important roles in various biological functions such as metabolism, immunity and development (Bartel, 2004;Nigam, 2014;Stepicheva and Song, 2015;Pan et al., 2021). Previous studies have recognized the critical role played by miRNAs during high temperature stress in aquatic animals (Li and Xu, 2018;Zheng et al., 2018;Ma et al., 2019). Studying gene and miRNA expression dynamics can illuminate the regulatory mechanisms of aquatic animals in response to environmental stresses. By combining RNA and miRNA analyses, a recent study identified a genetic advantage in the strong vigor group of Chinese mitten crab (Eriocheir sinensis) in terms of tolerance to nitrite stress. Meanwhile, several key genes involved in signal transduction, oxidative phosphorylation, and energy metabolism (such as HSP70, MTND2, ACT7, HSPG2, and Cam1) were detected. In addition, several miRNAs mediating calmodulin-related genes (miR-31, miR-215, and miR-142-x) were identified, which suggested that the calcium signaling pathway is involved in the response to nitrite stress in crabs (Fu et al., 2021). Qiang et al. (2017) constructed RNA and miRNA expression profiles of genetically improved tilapia farmed under high temperature stress and identified 28 differentially expressed (DE) miRNAs and 744 DE mRNAs as well as 64 negative miRNA-mRNA interactions. However, there has been no detailed investigation of miRNAs in response to heat stress in echinoderms, and currently, there is limited information on the effects of and the regulatory pathways involved in heat stress in S. intermedius.
This study is the first report of mRNA-miRNA interactions in S. intermedius under high temperature stress. We analyzed the physiological effects of heat shock on red and white tube-footed S. intermedius, detected and compared the dynamic changes of mRNAs and miRNAs in red and white tube-footed S. intermedius populations under high temperature stress, and provided clues to mRNA-miRNA interactions. These data help to understand the regulatory mechanism of S. intermedius on high temperature stress, provide information for screening molecular markers related to heat resistance and cultivating resistant sea urchins, and lay a foundation for healthy sea urchin farming. Before conducting the experiments, the sea urchins were acclimatized in two 0.5 m³ temperature-independent recirculating pools (water temperature maintained at 15 ± 0.5°C, salinity 31 ± 0.5‰, pH 8.2 ± 0.2) for 1 week. During the acclimation period, the sea urchins were fed fresh kelp daily, and the water was changed by 1/2, continuously aerated. After 48 h of fasting, the red and white tube-footed sea urchin population was randomly divided into 4 groups, namely, red tube-footed sea urchin control group (NR), white tube-footed sea urchin control group (NW), red tube-footed sea urchin high temperature group (HR), and white tube-footed sea urchin high temperature group (HW). The experimental design is shown in Figure 1.

Experimental Animals and Treatment
The control groups were the same as the acclimation groups, 15°C, and the experimental groups were 25°C (warmed at a rate of 2°C per hour), which was the lethal temperature for intermediate spherical sea urchins in our previous study. After reaching the temperature and maintaining it for 7 d, all sea urchins in each group were counted and individually weighed. Then, coelomic fluid of 15 randomly selected sea urchins from each group was extracted using 1 ml sterilized syringes. The upper fluid was obtained from the coelomic fluid by centrifugation (3,000 r/ min, 4°C, 10 min) for determining the activities of antioxidant enzymes. The lower precipitates were used for RNA extraction. All samples were frozen in liquid nitrogen and stored at −80°C.

Antioxidant Enzyme Activity Assays
Three sea urchin coelomic fluids were randomly selected from each group to detect antioxidant enzyme activities, including superoxide dismutase (SOD), catalase (CAT), peroxidase (POD) and glutathione peroxidase (GSH-Px), according to the instructions in the kit. All antioxidant enzyme test kits were purchased from the Nanjing Jiancheng Biological Engineering Institute (Nanjing, China).

RNA Library Construction and Sequencing
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 . High purity 1.5 μl of total RNA from each sample was interrupted, then was reverse-transcribed to construct cDNA libraries according to the instructions of the Illumina TruSeq RNA Sample Preparation Kit (Illumina, United States). Finally, paired-end sequencing was performed on an Illumina Novaseq ™ 6,000 (LC Sciences, United States). 12 mRNA libraries of sea urchins (three replicates of NR, NW, HR and HW) were built.

Mapping and Transcriptome Assembly
To explore the changes in gene expression profiles of different colored tube-footed sea urchin populations under hightemperature stress, we obtained high-quality clean data by removing missing bases and low-quality bases (base quality ≤10) (Martin, 2011), mapped it to the reference genome (not yet published) using hisat2, and assembled the transcriptomes of red and white tube-footed sea urchin coelomocyte as well as their coelomocyte under high-temperature stress using Trinity 2.4.0 (Grabherr et al., 2011). After the final transcriptome was generated, StringTie and ballgown (http://www.bioconductor. org/packages/release/bioc/html/ballgown.html) were used to estimate the expression levels of all transcripts and perform expression level for mRNAs by calculating FPKM (FPKM = [total_exon_fragments/mapped_reads (millions) × exon_length (kB)], command line:~stringtie -e -B -p 4 -G merged. gtf -o samples. gtf samples. bam).

Identification of Differentially Expressed Genes
Identification of DEGs was performed with edgeR software. The false discovery rate (FDR) was obtained by correcting the p-value using the Benjamini Hochberg's method. However, since very few genes were screened with FDR in this study, in order to expand the number of differential genes and find more differentially expressed genes, p ≤ 0.05 and |log 2 FoldChange| >1 were considered as the threshold to determine DEGs. All the DEGs were annotated according to GO (GOseq 1.10.0 method) and KEGG (KOBAS 2.0.12 method) (Mao et al., 2005;Young et al., 2010). p < 0.05 was set as a filtering threshold to detect significantly enriched GO entries and KEGG pathways associated with high temperature stress.

Small RNA Library Construction, Sequencing, and Annotation
The small RNA (sRNA) libraries were constructed with TruSeq small RNA sample prep kits (Illumina, San Diego, United States) and sequenced using Illumina HiSeq2500 (LC Sciences, Houston, Texas, United States). 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
To identify DE miRNAs between treatment group (HR, HW) and control group (NR, NW), the expression levels of miRNAs in the whole library were normalized, and then DE miRNAs were screened using Student's T-test. p < 0.05 was considered significant. 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
All possible positively and negatively correlated miRNA-mRNA pairs were predicted using ACGT101-CORR1.1. Based on the comprehensive analysis of DEG and DE miRNAs, we screened the negatively associated DE miRNA-mRNA pairs and constructed interaction networks using Cytoscape software (http://www.cytoscape.org/).

Validation of DE miRNAs and mRNAs Using qRT-PCR
To validate the results of RNA-seq and miRNA-seq, 10 DEGs and 10 DE miRNAs were randomly selected for real-time fuorescence quantitative PCR (qRT-PCR). Total RNA extraction and reverse transcription were performed by referring to Han et al. (2019). 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 Supplementary Table S1. The qRT-PCR was performed in a 20 μl reaction sample containing 2 μl of cDNA (50 ng/μl), 10 μl of FastStart Essential DNA Green Master, 6 μl of ddH 2 O water, and 1 μl (10 μM) 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 levels was determined using 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, United States). The data were expressed as the mean ± standard error (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.

Activities of Antioxidant Enzymes
Under high temperature stress, catalase (CAT) expression in the coelomic fluid of S. intermedius showed almost no change, whereas superoxide dismutase (SOD) and peroxidase (POD) expressions both increased significantly (p < 0.05), with greater changes in SOD in the white tube-footed S. intermedius and greater changes in POD in the red tubefooted S. intermedius. There was no significant change in glutathione peroxidase (GSH-PX) activity in the coelomic fluid of red tube-footed S. intermedius (p > 0.05), while the Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 876308 GSH-PX activity in the coelomic fluid of white tube-footed S. intermedius showed an increasing trend ( Figure 2).
The mapping results showed that 45.79-70.20% of clean reads were mapped to the reference genome sequence ( Supplementary Table S3), of which, on average, 80.20, 5.25, and 14.55% of the reads were mapped to the exon, intron and intergenic regions, respectively (Supplementary Table S4).

Identification and Functional Annotation of the Differentially Expressed Genes
On the basis of the sequencing data of red and white tube-footed S. intermedius in the normal group at 15°C (NR, NW) and under high temperature stress at 25°C (HR, HW), gene expression comparisons between control and treatment groups were undertaken.
In the HW vs. NW group comparison, a total of 407 DEGs (up-regulated: 178, down-regulated: 229) were detected, of which 296 DEGs were annotated in 1059 GO terms. The GO entries that are consistent with the HR vs. NR group were "microtubule-based process", "protein folding" and "transport"; "extracellular space", "extracellular exosome", "endoplasmic reticulum"; "calcium ion binding", "structural constituent of cytoskeleton" and so on. However, the difference was that "endoplasmic reticulum stress", "negative regulation of transcription" and "DNAtemplated" were found among BP and "RNA binding", "nuclide acid binding" and other structure-related terms were found in MF, while the DEGs were not annotated to the cell membrane entries in CC ( Figure 4B). KEGG enrichment analysis showed 137 DEGs were annotation to 108 pathways, and 17 were significantly enriched (p < 0.05). Different from the HR vs. NR group, more basal metabolism-related pathways were annotated, such as "one carbon pool by folate (ko00670)", "glyoxylate and dicarboxylate metabolism (ko00630)", and "linoleic acid metabolism (ko00591)" ( Figure 5B).
In the HR vs. HW group comparison, we detected 367 DEGs (up-regulated: 205, down-regulated: 162) of which there were 271 DEGs annotated to 851 GO terms. The most significantly enriched subcategories in BP, CC and MF were "microtubulebased process", "integral component of membrane" and "calcium ion binding", respectively ( Figure 4C). KEGG analysis showed a total of 118 DEGs were annotated in 83 pathways. The significant pathways containing the most DE genes were "phagosome (ko04145)", "ribosome (ko03010)", "galactose metabolism (ko00052)", "metabolism of xenobiotics (ko00980)", "cholinergic synapse (ko04725)" and so on. The top 20 enriched pathways for each comparison are shown in Figure 5C.

Identification and Functional Annotation of the Target Genes of the DE miRNAs
To explore DE miRNAs in the different comparison groups and further understand the regulatory functions of the miRNAs, we screened DE miRNAs and performed GO and KEGG enrichment analysis on the predicted target genes. A total of 2,543, 1,190, 1744, 984, 3,614, and 5,345 target genes were predicted for the six, two, three, one, five, and ten DE miRNAs identified from the six comparison groups, respectively (Supplementary Table S8), and GO analysis was used to identify enriched functional groups (p < 0.05). The top 10 significantly enriched GO terms (with lowest p-values) of the three categories are shown in Figure 7. Among the three groups, "regulation of transcription, DNAtemplated", "oxidation-reduction process", "cytoplasm", "nucleus", "ATP binding" and "protein binding" were the entries with the highest significant enrichment of target genes in BP, CC and MF.
KEGG pathway annotation showed 604, 316, 196, 316, 629 and 1,117 predicted target genes were annotated into 193, 161, 128, 161, 200, and 233 pathways, respectively (p < 0.05). Among them, more target genes were enriched in the classes of cellular processes, environmental information processing and metabolism. KEGG pathway analysis based on all of the predicted target genes revealed that DE miRNAs were involved in several pathways relevant to sea urchins exposed to high temperature, including protein processing in the endoplasmic reticulum (ko04144) and ubiquitin mediated proteolysis (ko04120), which are related to protein synthesis, and amino sugar and nucleotide sugar metabolism Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 876308 11 (ko00520) and alanine, aspartate and glutamate metabolism (ko00250), which are involved in basic metabolism. Furthermore, several immune-related pathways have also been discovered, including phagosome (ko04145), endocytosis (ko04144), autophagy (ko04140), mTOR signaling (ko04150) and FoxO signaling (ko04086) (Figure 8).

Prediction of Interactions Between DE miRNAs and DEGs
Combining the results for the DE miRNAs target genes and DEGs, we found that 5 DE miRNAs may have targeted regulatory relationships with 28 DEGs in the HR vs. NR group. Among these, 35 "DE miRNA-DEG" pairs exhibited a negative expression correlation, whereas 17 "DE miRNA-DEG" pairs showed a positive expression correlation. In the HW vs. NW group, we found that 1 DE miRNA may have targeted regulatory relationships with 20 DEGs. Among these, a negative expression correlation was found for 5 "DE miRNA-DEG" combinations, while a positive expression correlation was found for 15 "DE miRNA-DEG" pairs. In the HR vs. HW group, 3 DE miRNAs may have targeted regulatory connections with 27 DEGs. There was a negative expression correlation for 19 "DE miRNA-DEG" pairs, whereas there was a positive expression correlation for 8 "DE miRNA-DEG" pairs. In the NR vs. NW group, 1 DE miRNAs may have targeted regulatory connections with 10 DEGs. Among these, a negative expression correlation was found for 2 "DE miRNA-DEG" combinations, while a positive expression correlation was found for 8 "DE miRNA-DEG" pairs. In the HR vs. NW group, 5 DE miRNAs may have targeted regulatory connections with 16 DEGs. There was a negative expression correlation for 17 "DE miRNA-DEG" pairs, whereas there was a positive expression correlation for 18 "DE miRNA-DEG" pairs. In the HW vs. NR group, 10 DE miRNAs may have targeted regulatory connections with 107 DEGs. There was a negative expression correlation for 94 "DE miRNA-DEG" pairs, whereas there was a positive expression correlation for 70 "DE miRNA-DEG" pairs. This study focused on the response of sea urchins to high temperature, so we constructed mRNA-miRNA network maps for HR vs. NR, HW vs. NW, HR vs. HW groups. (Figure 9 and Supplementary Table S9). Since transcriptional attenuation of target genes is a conserved mechanism of miRNA regulation, we only focused on DEGs and DE miRNAs with negative regulatory relationships. For the "protein processing in endoplasmic reticulum (ko04141)" and "phagosome (ko04145)" pathways, we used Cytoscape software to map the network of DE miRNAs and targeted DEGs of these two pathways and found that miR-92b-3p and PC-5p-7420 were the most connected miRNAs, and lva-miR-92b-3p_R+2 and pmi-miR-92c-3p_R+1_1ss10TC co-regulate CALR (Table 2; Figure 10).

Validation of DEGs and DE miRNAs by qRT-PCR
qRT-PCR analysis was performed to validate the results of the differential gene expressions obtained from the RNA-Seq data. A total of 20 DEGs (such as HSP70, DnaJ11, and CALR) and 9 DE miRNAs were selected for the qRT-PCR measurements. The results of qRT-PCR analysis coincided with the results generated from RNA-Seq and showed similar FIGURE 9 | The network of DE miRNA and their target DEGs was constructed using Cytoscape 3.3.0. Yellow hexagon is miRNA, gray hexagon is target DEGs, red represents up-regulation, blue represents down-regulation.
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 876308 expression trends to those of the RNA-Seq data (Figures 11,  12). In addition, qRT-PCR results showed that the expression of HSP70, DnaJ11, GRP78, HYAL, and HYOU1 was significantly enhanced under high temperature stress, and the expression was higher in the red tube-footed sea urchin than in the white tube-footed sea urchin.

DISCUSSION
Temperature is the most important environmental element impacting organisms, and it can alter physiological and metabolic processes directly or indirectly. Mariculture organisms are primarily found in environmentally fragile coastal zones or estuaries, which are particularly vulnerable to environmental change (Dong et al., 2015;Ma et al., 2021). The S. intermedius is highly sensitive to environmental changes as a cold-water species, and significant temperature changes are one of the leading causes of mass mortality of sea urchins in China during the summer. Published biological and ecological data on the morphological forms of S. intermedius indicated that milkwhite color of spines (G) prevails in deep-water settlements (15-25 m) and the "usual" (U) form that mostly occurs in shallow-water settlements (5-10 m) (Balakirev et al., 2008;Balakirev et al., 2016). The milk-white color form preferring deep water conditions indeed could be less tolerant to hightemperature stress because the temperature fluctuations are much more significant in shallow-water than in deep-water environments (especially in the summer time). This study was conducted to investigate the changes in antioxidant enzyme activity under high temperature stress and to explore the candidate genes and regulatory pathways of S. intermedius in response to high temperature stress using RNA-Seq technology in order to better understand the molecular mechanisms invoked by sea urchins in response to high temperature stress. In response to environmental stress, organisms usually activate humoral immunity and alter their antioxidant enzyme activity to buffer the effects of temperature changes on the organism (Lu et al., 2016;Forgati et al., 2017). The antioxidant system, which primarily acts through antioxidant enzymes (such as SOD and CAT) and non-enzymatic antioxidants, protects organisms from oxidative damage (such as GSH). SOD is the first stage in the antioxidant process (Winston, 1991;Grundy and Storey, 1998), and its activity refers to the ability of an organism to scavenge free radicals, which is a key sign of its health (Luo et al., 2017). In this study, we evaluated the enzyme activities after 7 days of stress at 25°C. Compared with the control group at 15°C, the SOD, POD and GSH-Px activities of red and white tubefooted sea urchins were elevated, and the SOD activity of red tube-footed sea urchins was higher than that of white tube-footed sea urchins and the changes in their POD activity were significantly different (p < 0.05). This finding is consistent with the study of Ding et al. (2015), who investigated the effect of temperature on the activity of immune-related enzymes in S. intermedius with different colored tube feet, and our results confirmed the argument that the SOD activity of red tube-footed sea urchins responded to high temperature earlier than that of white tube-footed sea urchins. Meanwhile, we found no significant difference in CAT activity of the sea urchins under high temperature stress, suggesting that sea urchins use SOD more to regulate their ability to scavenge free radicals, or perhaps that enzyme activity is chronological and the magnitude of enzyme activity is related to sampling time. Therefore, we suggested that at the physiological-ecological level, red tubefooted sea urchins are more sensitive to high temperatures and have the potential to tolerate high temperatures.
In order to investigate the molecular mechanisms of sea urchins with varied tube foot colors in response to high temperature stress, gene expression profiles and miRNA profiles of S. intermedius with different tube foot colors were constructed in this study. Compared to previous studies (Ding et al., 2017;Zhao C et al., 2018), this study used our own sequenced genome of S. intermedius to construct a higher quality of data with the reference transcriptome and screened for more DEGs. We found that under normal rearing, the immune-related pathways such as "phagosome" and "lysosome" were significantly enriched in the red tube-footed S. intermedius. We found the number of DEGs in red tube-footed sea urchins was less than that in white tube-footed sea urchins under high temperature stress compared with the control groups.
FIGURE 10 | The integrating biomolecular interaction network of target DEGs predicted to be regulated by DE miRNA in the "protein processing in endoplasmic reticulum (ko04141)" and "phagosome (ko04145)" pathways. Red diamond is KEGG pathway, yellow hexagon is miRNA, gray hexagon is target DEGs, red represents up-regulation, blue represents down-regulation.
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 876308 FIGURE 11 | Analysis of selected 20 differentially expressed genes (DEGs) in red and white tube-footed S. intermedius under the high temperature stress by sequencing and qRT-PCR. Each vertical bar represents the Mean ± SD (n = 3), 18s rRNA were used as a reference gene. *Significant differences at p < 0.05 vs. control (NR and NW). **Highly significant differences at p < 0.01 vs. control (NR and NW).
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 876308 16 This finding is in accordance with earlier research (Shi et al., 2020). 20 DEGs and 9 DE miRNAs were screened for qRT-PCR analysis to validate the RNA-Seq results, and the results showed that the qRT-PCR expression trends were consistent with the RNA-Seq data, validating the reproducibility of our RNA-Seq research.
The "protein processing in endoplasmic reticulum" pathway, which alters protein folding in the endoplasmic reticulum, was the pathway most affected by high temperature stress , according to the KEGG enrichment analysis of DEGs in the HR vs. NR and HW vs. NW groups. The expressions of HSP70, GRP78 and PDIA6 were considerably elevated under high temperature stress, showing that heat stress causes protein misfolding and triggers the stress response. This result is consistent with the results of Yang et al. (2017), who studied the heat stress response mechanism in oysters (Crassostrea gigas). Heat shock proteins are chaperone proteins that promote the degradation of abnormal proteins or reactivate stress-damaged proteins by promoting the proper refolding of denatured proteins (Morimoto et al., 1992;Zhao Z. W et al., 2018). The increased expression of HSP70 would improve the thermal tolerance of cells under high temperature induction, limiting the production of oxygen radicals caused by heat stress and thereby reducing the rate of heat stress-related cell death (Zhou et al., 2004). Previous research has found that S. intermedius HSP70 responds to rapid temperature changes in a variety of tissues (coelomocyte, tube feet, intestine, gonads) with differential expression and that HSP70 expression increases and then decreases with stress time (Bai et al., 2015). GRP78 (glucose-regulated protein 78 kDa) is a member of the HSP70 family, and it is a key indicator of endoplasmic reticulum stress (Lee, 2005). This reveals that these DEGs may help misfolded proteins refold into the right molecular shape. The expressions of essential unfolded protein response (UPR) components including calreticulin (CALR) and protein disulfide isomerase (PDI) were also shown to be significantly up-regulated. CALR is a Ca 2+binding protein that regulates intracellular Ca 2+ dynamic homeostasis and signaling activities such as the UPR and related protein degradation pathways, as well as regulating the immunological response, shielding cells from stress and regulating the cell cycle (Eggleton et al., 2016). PDI promotes protein folding by acting as a chaperone molecule (Schröder and Kaufman, 2005). Ding et al. discovered 18 distinct immunerelated pathways in the transcriptome of red tube-footed sea urchins when compared to white tube-footed sea urchins, implying that red tube-footed sea urchins have a stronger immune response than white tube-footed sea urchins (Ding et al., 2017). We also discovered 13 potentially specific KEGG pathways of red tube-footed sea urchins under high temperature stress, including "phagosome", "ribosome" and "cholinergic synapse", as well as other immune-related and basal metabolic pathways such as "gluconeogenesis". The immunological FIGURE 12 | qRT-PCR verification of select 9 DEMs. Each vertical bar represents the Mean ± SD (n = 3), 18s rRNA were used as a reference gene. *Significant differences at p < 0.05 vs. control (NR or NW). **Highly significant differences at p < 0.01 vs. control (NR or NW).
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 876308 response and basal metabolism of red tube-footed sea urchins are thought to be boosted in high-temperature environments to protect the organism from stress damage. It can thus be suggested that red tube-footed sea urchins are more sensitive than white tube-footed sea urchins to respond to heat stress. Gene expression is thought to be influenced by a variety of factors. Recently, researchers have investigated changed gene expression at the post-transcriptional stage in a variety of biological systems using miRNA-mRNA regulatory networks (Zhou et al., 2019;Wang et al., 2020). The miRNA profiles of the sea urchin S. intermedius under temperature stress are presented in this study. The experimental results demonstrated that miRNAs are mostly 22 nt in length, which is consistent with earlier research (Zhan et al., 2018). The 11 identified DE-miRNAs were down-regulated and involved in multiple key biological processes related to biosynthesis, metabolism, immunity and signaling transduction.
Since transcriptional attenuation of target genes is a conserved mechanism of miRNA regulation, we only focused on "DEG-DE miRNA" pairs with opposite relative expression trends in this study. After screening, we found that these "DEG-DE miRNA" pairs may explain the superiority of energy metabolism and immunity of red tube-footed sea urchins in response to high temperature stress. It has been shown that miR-184 is an important regulator of stem cell proliferation and growth , and that its overexpression can lead to apoptosis and its suppression can increase cell numbers (Foley et al., 2010). Previous studies have shown that miR-184 overexpression inhibits autophagy and exacerbates oxidative damage, and it can also negatively regulate Wnt signaling in vivo and in vitro . MiR-184 could inhibit protein expression in human trabecular meshwork cell cytotoxicity, apoptosis and the extracellular matrix via targeting HIF-1α in vivo, and it can also exhibit angiostatic properties through regulating signaling pathways including Akt, TNF-α and VEGF (Takahashi et al., 2015). The target gene of miR-184 was discovered to be hypoxia up-regulated protein 1 (HYOU1) in this study, and the two had a negative regulatory connection. When cells are stressed by the external environment (hypoxia, ischemia and high temperature etc.), HYOU1, as an endoplasmic reticulum molecular chaperone, protects and maintains cell viability (Park et al., 2017). Therefore, we hypothesized that miR-184 regulates the expression of genes involved in cell proliferation and autophagy in red tube-footed sea urchins in response to heat stress. We also found that red tube-footed sea urchins displayed downregulation of miR-92a, miR-92c and miR-92b under high temperature stress. MiR-92a belongs to the miR-17-92 cluster, which has six mature miRNAs (Tsuchida et al., 2011). Members of this cluster are involved in the development of vascular endothelial cells and regulate numerous essential biological processes such as cell proliferation and death (Tamatani et al., 2001;Olive et al., 2009;Olive et al., 2010;Kaluza et al., 2013;Tan et al., 2014). Additionally, miR-92a promotes cell proliferation by inhibiting an isoform of the cell cycle regulator p63 (Manni et al., 2009). Yang et al. also showed that miR-92d can regulate the complement system in the amphioxus Branchiostoma belcheri to modulate the acute immunological response to bacterial infections (Yang et al., 2013).  is involved in apoptotic signaling pathway regulation perhaps via targeting Aj14-3-3ζ in sea cucumbers (Lv et al., 2017). MiR-92c and miR-92b were discovered to influence CRT expression in our study, and more research will be done to confirm this. We also discovered miR-124 by comparing the miRNA expression profiles of red and white tube-footed sea urchins at high temperatures. In the mouse brain, miR-124 promotes neuronal development by triggering particular alternative pre-mRNA splicing (Makeyev et al., 2007). Non-neuronal transcripts are increased in primary cortical neurons when miR-124 is inhibited (Conaco et al., 2006). We hypothesized that decreased miR-124 expression led to a stimulation of neural differentiation in red tube-footed sea urchins under high temperature stress, making them more temperature sensitive and thus faster to respond to high temperature stress. Two new miRNAs, PC-3p-45817 and PC-5p-7420, were identified by association analysis, in which the target genes of PC-2p-4581 were DNAJ11 and PDIA6, and the target genes of PC-5p-7420 were GANAB, SAR1B and TUBA. These genes play an important role in the UPR response to endoplasmic reticulum stress (Gotoh et al., 2004;Sha et al., 2012). Our findings support the idea that miRNAtarget gene interactions are complicated, and that a single miRNA can regulate multiple target genes simultaneously, whereas a single gene can be regulated by multiple miRNAs at the same time (Lan et al., 2016). The results of the qRT-PCR validation further imply that there may be meaningful targeted regulatory interactions between th candidate DE miRNAs and DEGs. As a consequence, in future research, we will need to confirm the targeted regulatory interactions of these putative "DEG-DE miRNA" pairs utilizing additional in vivo and in vitro experiments (e.g., dual luciferase reporter analysis, gain or loss of function analysis and Western blotting).

CONCLUSION
To research the molecular processes involved in the sea urchin response to high thermal stress, we investigated gene expression differences and miRNA profiles of red and white tube-footed sea urchins under high temperature stress conditions. According to previous studies and miRNA-mRNA integration analysis, protein folding of the endoplasmic reticulum and phagosomes are the main pathways for sea urchins to respond to high temperature stress, and red tube-footed sea urchins may be more sensitive to high temperature and respond more quickly to high temperature stress. The red tube-footed sea urchin may in fact be more sensitive to high temperatures than other sea urchins and respond faster to heat stress. Red tube-footed sea urchins and "DEG-DE miRNA" pairs involved in protein processing of the endoplasmic reticulum may be considered as germplasm resources and molecular markers for future selective breeding.

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 in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
JD and YC: conceptualization and resources. JD and LH: conceived and designed the experiment. ZQ, LH, PH, and WW: performed the experiment. LH, HW, and LW: data curation. LH, XZ, PL, and YL: analyzed the data. YW, CG, DY, and WZ: contributed reagents/materials/analysis tools. ZQ 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 National Key R&D Program of China (2018YFD0901601), High-level Talent Support Grant for Innovation in Dalian (2020RD03), and Liaoning Province Higher Education Innovation Team and Innovative Talent Support Program Project (LT2019003).