Differential expression of miRNAs in the body wall of the sea cucumber Apostichopus japonicus under heat stress

MicroRNAs, as one of the post-transcriptional regulation of genes, play an important role in the development process, cell differentiation and immune defense. The sea cucumber Apostichopus japonicus is an important cold-water species, known for its excellent nutritional and economic value, which usually encounters heat stress that affects its growth and leads to significant economic losses. However, there are few studies about the effect of miRNAs on heat stress in sea cucumbers. In this study, high-throughput sequencing was used to analyze miRNA expression in the body wall of sea cucumber between the control group (CS) and the heat stress group (HS). A total of 403 known miRNAs and 75 novel miRNAs were identified, of which 13 miRNAs were identified as significantly differentially expressed miRNAs (DEMs) in response to heat stress. A total of 16,563 target genes of DEMs were predicted, and 101 inversely correlated target genes that were potentially regulated by miRNAs in response to heat stress of sea cucumbers were obtained. Based on these results, miRNA-mRNA regulatory networks were constructed. The expression results of high-throughput sequencing were validated in nine DEMs and four differentially expressed genes (DEGs) by quantitative real-time polymerase chain reaction (RT-qPCR). Moreover, pathway enrichment of target genes suggested that several important regulatory pathways may play an important role in the heat stress process of sea cucumber, including ubiquitin-mediated proteolysis, notch single pathway and endocytosis. These results will provide basic data for future studies in miRNA regulation and molecular adaptive mechanisms of sea cucumbers under heat stress.

Introduction

The sea cucumber (Apostichopus japonicus) is widely distributed along the northwest Pacific coast at 35∼44°N and has been utilized as an important marine economic resource for its nutritional and medicinal value (Chang and Loscalzo, 2010; Wang et al., 2015). The total amount of sea cucumber cultivation in China has continuously increased and the cultivation area has been expanding since 2003. In 2020, the annual production of sea cucumber cultivation reached 200,000 t with an increase of 14.48% over the previous year, and the production value reached 50 billion US dollars (Ministry of Agriculture and Rural Affairs, 2021). The sea cucumber aquaculture industry has suffered great loss, with the global warming and the frequent occurrence of high temperatures in summer. Heat stress was reported to impose many negative effects on the growth, metabolism, and reproduction of aquatic farmed animals. Therefore, understanding the molecular mechanism of miRNA under heat stress would help to reduce the damage caused by heat stress on sea cucumbers and breed new varieties with high temperature resistance.

As a typical temperate species, the sea cucumber was greatly affected by water temperature, not only on its physiological activities, but also on feeding and metabolism (Zhao et al., 2015; Li et al., 2016a). The high temperature was shown to reduce the activity and food intake of sea cucumber, disorder the metabolism of free radicals in the body and even lead to disease and death (Wang et al., 2011; Li et al., 2012). It was showed that the total number of coelomic cells, phagocytosis rate and lysozyme activity in the coelomic fluid of A. japonicus was decreased to a certain extent after entering summer dormancy (Shao et al., 2016). The research on the effects of temperature changes on the physiology, growth and gene expression of sea cucumbers has gotten a lot of attention. Zhang et al. (2015) found that increasing temperature can affect the DNA methylation level of A. japonicus, and then change the expression of the heat shock protein gene. The methylation sites of intestinal, respiratory tree and gonadal tissues in sea cucumber were mainly distributed in the gene functional regions, and the methylation sites were mainly of CG type (Li et al., 2018). Under heat stress at 26°C, the methylation level of intestinal genome in sea cucumber increased, while at 32°C, the methylation level decreased (Wen et al., 2021). At present, the research on the regulation mechanism of sea cucumber under heat stress focuses on the epigenetic and transcriptional regulation, but there are few studies on the post-transcriptional regulation.

MicroRNAs (miRNAs) are endogenous noncoding small RNAs with 22 nt, which regulate post-transcriptional gene expression by binding to the 3′ untranslated region (UTR) of target genes (Bartel, 2004; Pedersen et al., 2007). Numerous studies have reported that miRNAs can provide reversible gene silencing mechanisms during animal aestivation and hibernation, thus miRNAs can be involved in cellular adaptation to specific demands under stressful conditions (Jones-Rhoades and Bartel, 2004). It has been reported that miRNA played a key role in the response to heat stress of plants and animals (Jones-Rhoades and Bartel, 2004). In rainbow trout, some miRNAs, including ssa-miR-301a-3p, ssa-miR-30a-5p and ssa-miR-30a-5p, can regulate the key changes of cells under high temperature stress (Ma et al., 2019). In Litopenaeus vannamei, 41 differentially expressed miRNAs related to heat stress were selected, such as lva-miR-92b, lva-miR-317 and lva-miR-184 (Boonchuen et al., 2020). Recently, Li et al. (2016b) and Zhou et al. (2018) found that some miRNAs including miR-184 and miR-2004 played a key role in the coelomic fluid of sea cucumber in response to heat stress. However, it was largely unknown for the miRNA regulation in the body wall of diseased sea cucumber in response to heat stress.

The purpose of this study was to identify known miRNAs and novel miRNAs from the body wall of sea cucumber high-throughput sequencing analysis, as well as to investigate the regulation mechanism between miRNAs and their target genes in sea cucumber after heat treatment. Our findings will be useful in further studies of sea cucumber biomarkers and will provide basic data for future studies on miRNA regulation and sea cucumber molecular adaptive mechanisms under heat stress.

Materials and methods

Experimental design and tissue collection

Sea cucumbers (body weight 50.0 ± 2.0 g) of appr. 1 year old were collected from Qingdao Ruizi Group Co., Ltd., in Shandong Province, China. The sea cucumbers were transported to our lab in Qingdao and acclimated in aerated indoor tanks for 5 days. During the experiment, the temperature and salinity of sea water were around 15°C and 30‰, respectively. The sea cucumbers were fed with a regular compound feed, and one third of seawater was changed daily.

When acclimation finished, a total of 30 sea cucumbers were selected and randomly divided into two groups. For the control group (CS), the culture temperature was continued to maintain at 15°C. According to previous studies of our lab, the median lethal temperature (LT₅₀) to sea cucumber was 32°C (Zhang et al., 2022). Therefore, the temperature for the heat stress group (HS) was set to 32°C in the present study. Heaters were used to give thermal stress to the sea cucumbers. The temperature system was continuously increased at a heating rate of 2°C/24 h. The moment when the water temperature reached 32°C was regarded as the initial time. In the subsequent experiment, the water temperature was maintained at 32°C ± 0.5°C. On the third day of the experiment, three sea cucumbers in the CS group and three skin ulceration syndrome individuals in the HS group were randomly selected. The body wall of six sea cucumbers were immediately sampled and frozen in liquid nitrogen, since the body wall was the target organ of skin ulceration syndrome. All the samples were stored at −80°C for miRNA sequencing analysis.

RNA extraction and processing

Total RNA was extracted from each sample using the Animal Tissue RNA Purification Kit (LC Sciences, Houston, TX, United States) according to the manufacturer’s instructions. The quality of total RNA was checked in Bioanalyzer 2100 (Agilent, Santa Clara, CA, United States) with RNA integrity number >7.0.

Small RNA library was constructed according to the protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, United States). RNAs in the 16–30 nt size range were purified from a 15% polyacrylamide gel, then sequentially ligated to 5′ and 3′ adapters. Reverse transcription was subsequently performed by polymerase chain reaction (PCR) amplification. The purified PCR products were sequenced by Illumina Hiseq2500 (LC-BIO, Hangzhou, China).

Sequence data analysis

Raw reads were subjected to an in-house program, ACGT101-miR (LC Sciences, Houston, TX, United States) to remove adapter dimers, junk, low complexity, and clean reads ranging from 16 to 30 nt in length were obtained. Sequence matching non-coding RNAs, including rRNA, tRNA, small nuclear RNA (snRNA), and small nucleolus RNA (snoRNA), were mapped to Repbase database (http://www.girinst.org) and Rfam database (http://rfam.xfam.org). The unique small RNA sequences were aligned against the miRBase V22.0 (http://www.mirbase.org/blog) to select the known miRNAs and novel miRNAs and at most one mismatch inside of the sequence was allowed in the alignment. The Deuterostoma species in the miRBase were used as host species and the priority rank of these species in sequence alignment was listed in Supplementary Table S1. Then, the sequence was mapped to the sea cucumber genome (assembled by our lab, unpublished) by Bowtie 1 (Berthelot et al., 2014) to determine the genomic locations of known miRNAs and the flank sequences of novel miRNAs. The hairpin structures of the novel miRNAs were predicated from the flank 80 nt sequences using miRDeep2 (Friedlnder et al., 2012). Principal component analysis (PCA) was performed on the valid data using the vegan R package.

TPM (millions of transcripts) normalization was used for miRNA expression analysis. Differential expression of miRNAs was selected using DEGseq (http://www.bioconductor.org/packages/release/bioc/html/DEG-seq.html). Fold-change and p-values were calculated from the normalized expression referring to the methods in published studies (Song et al., 2017). |log₂(Foldchange)| > 1 and p-value lower than 0.05 was considered as significantly different expression miRNAs (DEMs).

Target gene prediction and function analysis

MiRanda (http://www.microrna.org/) and TargetScan (http://www.targetscan.org/vert_80/) were used to predict the target genes of known miRNAs and novel miRNAs. Target genes with a context score percentile <50 were removed from the TargetScan, and target genes with maximum free energy (Max Energy) >−10 were removed from the miRanda. Finally, the intersection of these two websites was taken as the final target genes of different expression miRNA.

Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) analyses were carried out to further understand genes biological functions. Target genes of differential expressed miRNAs were mapped to GO terms in the database (http://www.geneontology.org/) and KEGG database (http://www.genome.jp/kegg/) respectively. The GO terms and KEGG pathways with p-value < 0.05 through hypergeometric test were defined as significantly enriched terms.

Correlation analysis of miRNA and differentially expressed genes

In previous mRNA-seq studies, we have identified differentially expressed genes (DEGs) under heat stress using the same samples as the present study (Li, 2021). Therefore, correlation analysis of miRNA-mRNA were carried out in order to identify key miRNA-target pairs. p < 0.05 and fold change >2 were set as the threshold for screening miRNA-target pairs. Only inversely correlated miRNA-mRNA were identified. The identified miRNA-mRNA regulatory network was constructed by Cytoscape 2.8.3 software (Shannon et al., 2003).

Quantitative real-time PCR validation of miRNA and mRNA expression

Quantitative real-time PCR (RT-qPCR) was used to verify the expression level of nine DEMs and four DEGs that were selected from miRNA-mRNA regulatory network based on their potential functional importance. The remaining RNA samples from the small RNA-Seq library construction were used for RT-qPCR amplifications. β-Actin and U6 were used as the internal control gene for DEGs and DEMs (Table 1). Total RNA after genome DNA removing was reverse-transcribed to cDNA using Evo M-MLV RT Kit (TaKaRa, Japan) following the manufacturer’s instructions. The reaction was set for 25 min at 37°C and 5 s at 85°C. The cDNA was amplified using SYBR Green Premix Pro Taq HS qPCR Kit (TaKaRa, Japan) with miRNA and mRNA specific forward primers and universal reverse primers. The PCR amplification parameters were as follows: 95°C for 2 min, followed by 45 cycles of 95°C for 15 s, 60°C for 15 s, and 70°C for 25 s. Amplification specificity was verified by melting curve analysis. The relative expression levels of target gene transcripts were calculated according to the comparative cycle threshold (Ct) method (2^−ΔΔCT), and GraphPad 6.0 was used for statistical analysis. All data are expressed as mean ± standard deviation (SD).

TABLE 1

TABLE 1. RT-qPCR primers used for validation of this study.

Results

Small RNA library construction

In total, six small RNA libraries (CS1, CS2, CS3, HS1, HS2, HS3) were constructed from the body wall of sea cucumbers in the present study. The raw data have been submitted to the National Center for Biotechnology Information under the accession number SRR18918495 ∼ SRR18918500. A total of 11,929,100 and 10,107,557 raw reads were generated from the CS and HS, respectively. After removing low-quality reads, including repeat reads, junk reads, mRNA and other noncoding RNAs, a total of 5,902,588 and 5,053,528 valid reads were obtained from the CS and HS, respectively (Table 2). PCA analysis was performed on the data of six groups, and it was found that HS group and CS group had great repeatability, respectively (Supplementary Figure S1). From the valid reads length statistics, the distribution showed that 22 nt was the most common type (Supplementary Figure S2). A total of 403 known miRNAs were identified in six libraries (Supplementary Table S1). Among these known miRNAs, 281 miRNAs were grouped into 127 families, of which miR-25 was the largest family with 16 members, followed by let 7 family with 10 members, and the miR-1 with eight members (Supplementary Table S3). Among the 127 families, 61 families contained only one number, and 33 families contained two numbers. The information of top 20 families were displayed in a histogram (Supplementary Figure S3). The sequences nonmatched to miRBase database were compared with the genome to predict the novel miRNAs by miRDeep2 and 75 novel candidate miRNAs were predicted (Supplementary Table S4). Some novel miRNAs with the same mature sequence but different precursors were considered to belong to a novel miRNA family. These novel miRNA candidates were named in the form of “PC plus number” (Supplementary Table S4). The miRNA precursor sequences plus a sequence of 100 nt per site were presented in Supplementary Table S5.

TABLE 2

TABLE 2. Statistics of small RNA sequences from the six libraries in the body wall of Apostichopus japonicas.

Differential expression analysis of miRNAs

Compared with the control group, 144 miRNAs were upregulated including miR-92a and miR-210, while 334 miRNAs were downregulated including miR-152 and miR-22. The t-test was used for analysis and p < 0.05 was regarded as the criterion to identify the differential expressed miRNAs (DEMs). The result showed that 13 DEMs were identified between CS and HS, and the numbers of upregulated and downregulated DEMs were five and eight, respectively (Figure 1; Table 3). These miRNAs may play an important role in regulating sea cucumbers in response to heat stress. Besides, there was a novel miRNA (Novel-6338) that included in these DEMs. Among these DEMs, the expression of miR-92a-3p and miR-210 showed more than two-fold up regulation and miR-10 showed more than five-fold up regulation. Particularly, the expression levels of miR-31, miR-2004-5p and miR-1357-3p were downregulated over four fold.

FIGURE 1

FIGURE 1. Cluster analysis heatmap of differentially expressed miRNAs.

TABLE 3

TABLE 3. Information of significantly differential expressed miRNAs in the body wall of Apostichopus japonicas under heat stress.

Target gene prediction and functional analysis

In order to understand the biological processes and molecular functions of miRNAs during heat stress, Targetscan and miRanda were used to predict the target genes of the 13 DEMs. TargetScan algorithm removes the target gene whose context score percentile is less than 50, and miRanda algorithm removes the target gene whose maximum energy is greater than −10. Finally, the intersection of these two softwares is taken as the final target gene of DEMs. For the 13 DEMs, a total of 16,563 target genes were predicted, including 5,153 target genes of five upregulated DEMs and 11,410 target genes of eight downregulated DEMs. The GO database annotation results showed that 11,001 target genes of 13 DEMs were annotated to 334 GO terms (p < 0.05), including 178 terms for biological processes, 65 terms for cellular components, and 91 terms for molecular functions. Of these GO terms, membrane, integral component of membrane and nucleus in cellular components are enriched with more genes, and we listed the top 10 terms with the most target genes number in each category (Figure 2). The result showed that 4,235 target genes were annotated into 25 KEGG pathways (p < 0.05), among which the endocytosis, ubiquitin mediated proteolysis, homologous recombination, axon guidance, notch signaling pathway interaction were the top five pathways with the most abundant target genes (Figure 3).

FIGURE 2

FIGURE 2. Top 10 enriched GO terms of the target genes of differentially expressed miRNAs.

FIGURE 3

FIGURE 3. KEGG pathway enrichment of the target genes of differentially expressed miRNAs.

Correlation analysis of miRNA-mRNA

To reveal the roles of miRNAs, the DEMs were selected for association analysis with the DEGs obtained by previously transcriptome sequencing. Finally, we identified 101 inversely correlated target genes that potentially regulated by miRNAs in response to heat stress of sea cucumbers (Supplementary Table S6). Of these 101 inversely correlated target genes, 12 inversely correlated target genes were down regulated by four up regulated DEMs, and 89 inversely correlated target genes were up regulated by eight down regulated DEMs. Among them, miR-92c combined with the most DEGs in the upregulated DEMs, and miR-2004-5p combined with the most DEGs in the downregulated DEMs (Figures 4A,B). The DEGs in the regulation network, such as WD repeat-containing protein 20 (WD20), glycoprotein 3-alpha-L-fucosyltransferase A-like (FUT4), E3 ubiquitin-protein ligase (TRAF), Serine/threonine protein kinase TAO3-like isoform (PAK4), interleukin -17D-like (IN17) and BAG family molecular chaperone regulator 3 (BAG3) involved in ubiquitin-mediated proteolysis and notch signaling pathway. These results suggested that these negative regulatory genes might be involved in the response of A. japonicus to heat stress.

FIGURE 4

FIGURE 4. Regulation networks of miRNA-mRNA of key genes of correlation analysis. (A) Down-up mode; (B) Up-down mode.

Validation of differentially expressed miRNAs and their differentially expressed genes expression by real-time PCR

To verify the reliability of the sequencing expression profiles, nine miRNAs and four negative genes were applied to RT-qPCR. A peak was detected in the melting curve during the experiment, indicating that all PCR products were specifically amplified. The results of RT-qPCR showed that the expression trend of miRNA and negative genes were consistent with the high-throughput sequencing results. The expression of miR-210 and miR-10, miR-92a-3p, miR-92a-5p and miR-92c were significantly upregulated (Figure 5A), and miR-2004-5p, miR-31, miR-193 and miR-4185 were significantly downregulated (Figure 5B). The expression of BAG3 and TARF7 were significantly upregulated, and WDR20 and FUT4 were significantly downregulated (Figure 5C).

FIGURE 5

FIGURE 5. Validation of differentially expressed miRNAs and target genes via RT-qPCR method; (A) Validation of upregulated differentially expressed miRNAs; (B) Validation of downregulated differentially expressed target genes; (C) Validation of negative DEGs. Values indicate the mean ± SD (n = 3).

Discussion

In the present study, we performed a high-throughput miRNA sequencing and miRNA-mRNA correlation analysis to study the post-transcription regulation mechanism of A. japonicus in response to high-temperature stress. We obtained 13 significantly differentially expressed miRNAs and a list of negative regulator genes of these differentially expressed miRNAs in the body wall of A. japonicus during heat stress. These miRNAs, negative target genes and miRNAs-mRNA regulation networks may play important roles in regulating heat stress, which broadens our understanding of the molecular regulation mechanism of A. japonicus in response to heat stress.

Among the 13 identified differentially expressed miRNAs, some miRNAs (such as miR-184, miR-193, miR-2004, miR-31, miR-210 and miR-10) have been reported to be involved in stress responses or immune processes previously. MiR-184 has been reported to be related to the infection of Procambarus clarkii by Spiroplasma eriocheiris (Ou et al., 2013). Xu et al. (2019) found that miR-184 of Sinopotamon henanense was significantly upregulated in response to oxidative stress induced by cadmium (Cd). MiR-193 has been reported to be related to Bos taurus infection (Shi et al., 2014). The expression level of miR-2004 was increased with the increase of A. japonicus complement AJC3 in the coelomic cells of A. japonicus in different periods after LPS stimulation (Zhong et al., 2015). The results of Zhou et al. (2018) showed that miR-2004 was regarded as a significantly expressed gene in the coelomic cells of A. japonicus under the infection of V. splendidus. MiR-31 has been proved to regulate the occurrence of inflammatory bowel disease by directly targeting the hypoxia-inducible factor (Olaru et al., 2015). Suzarez et al. (2010) confirmed that miR-31 could regulate the inflammatory response by negatively regulating the binding of e-selectin and neutrophils to endothelial cells. Kang et al. (2019) found that miR-31 was significantly upregulated in the liver after infecting grouper with V. alginolyticus, and the target genes were mainly concentrated in immune-related pathways. Therefore, we inferred that miRNA plays a complex regulatory role in various stress responses.

MiR-210 was known as a major hypoxia-inducible microRNA, which is evolutionarily conserved and ubiquitously expressed in the hypoxic cell and tissue types, serving versatile functions (Chang and Loscalzo, 2010). It was found that miR-210 can negatively regulate the production of LPS-induced pro-inflammatory cytokines through NF-κB1 in mouse macrophages (Qi et al., 2012). Yu et al. (2021) found that miR-210 can activate PI3K/AKT pathway by regulating PI3K and p-AKT protein, promoting the proliferation, and inhibit the apoptosis of dental pulp stem cells. Recently, Li et al. (2016b) found that miR-210 regulated the host defense of A. japonicus through TLRs (Toll-like receptors) pathway. In the present study, the significantly upregulated expression indicates that miR-210 also play role in sea cucumber’s body wall in response to heat stress. It was speculated that A. japonicus may reduce the expression of target genes by upregulating the expression level of miR-210 to cope with the harm of heat stress.

MiR-10 family members were highly conservative and closely associated with metabolizing response. MiR-10 can directly regulate the expression levels of FLT1 (vascular endothelial growth factor receptor 1) and sFLT1 (soluble vascular endothelial growth factor receptor 1), both of which are closely related to the growth of vascular endothelial cells (Livak and Schmittgen, 2013). Knock-out of miR-10 could lead to premature termination of the growth of blood vessels in the embryonic internodes of zebrafish larvae, while over-expression of miR-10 could accelerate the angiogenesis in zebrafish and the growth of human umbilical vein endothelial cells (Wienholds, 2005). Another study, found that miR-10 can regulate the growth of blood vessel endothelial cells by promoting new signal transmission (Naguibneva et al., 2006). Overexpression of miR-10 could promote virus-induced apoptosis, and researchers discovered that the expression of apoptosis protein caspase-3 increased in tandem with miR-10 expression (Zhang, 2014). Moreover, miR-10 was significantly upregulated in the process of infection by infecting porcine alveolar macrophages with PRRSV (Porcine Reproductive and Respiratory Syndrome Virus), and the transfection experiment demonstrated that the upregulated expression of miR-10 could significantly inhibit the replication of HP-PRRSV and N-PRRSV in porcine alveolar macrophages (Zhao et al., 2017). In the present study, we found that miR-10 was significantly upregulated 5-fold after heat stress. It indicates that miR-10 regulates the metabolism of A. japonicus during heat stress, but it is molecular function needs further to be studied.

In the process of A. japonicus in response to heat stress, ubiquitin-mediated proteolysis and notch signaling pathway were significantly enriched, indicating that they played important role in the process. Ubiquitin-mediated proteolysis is an important cellular immune pathway, which can control the basic life activities of cells by degrading key regulatory proteins and regulating the cellular stress response and immune response to pathogenic microorganisms (Chen et al., 2012). For example, when shrimp was infected by V. cholerae, ubiquitin proteins can mediate cellular immunity to cope with the infection (Li, 2020). The notch signaling pathway can regulate cell differentiation, proliferation and apoptosis, which is of great significance to cell growth and development. Abnormal expression of notch signaling can induce many kinds of cancers, such as breast cancer (Wang et al., 2017), lung cancer (Tian et al., 2017), and gastric cancer (Wang et al., 2018). We inferred that ubiquitin-mediated proteolysis and Notch signaling pathway could resist the effects of high temperature by regulating the expression of downstream immune genes in the process of heat stress.

Through the correlation analysis between the differentially expressed miRNAs and transcriptome data, it found that 12 DEMs can target 101 inversely correlated DEGs. TRAF7 has been partially studied in aquatic animals, and BAG3, WDR20 and FUT4 are mostly focused on wound tissue recovery and cancer treatment in medicine. BAG3 is a member of the BAG gene family, which plays an important role in regulating cell apoptosis, autophagy, movement and development, and mediating the adaptability of cells response to heat stress (Rosati et al., 2011). The results of GO enrichment analysis showed that BAG3 could inhibit cell apoptosis and promote tumor cell proliferation by binding to HSP70 (Manzerra and Brown, 1990). TRAF7 belongs to the TRAFs family, which involves a variety of biological functions including innate immunity, embryo development, stress response and inflammatory response (Reuss et al., 2013). In mammals, TRAF7 can regulate the signal activities by enhancing MEKK (mitogen-activated protein kinase), and then participate in the immune response of the body (Xu et al., 2004). Cynoglossus semilaevis has also been found to downregulate TRAF7 expression in all tissues infected with V. harveyi, particularly the liver (Wei et al., 2018). Studies have shown that WDR20 has a positive or negative correlation with the occurrence of lymphoma, prostate cancer and other diseases (Ohashi et al., 2015). At the same time, experiments have also proved that miR-3188 can target and inhibit the expression of WDR20 (Zhang et al., 2018). Screening of key genes of A. japonicus in response to heat stress provides basic data for elucidating its regulation mechanism. Further research will focus on the verification of the regulatory relationship between DEMs and DEGs.

Conclusion

In this study, we identified the expression profiles of miRNAs under heat stress in A. japonicus by using small RNA-seq. We highlighted 13 DEMs compared with the control group, which were involved of immunity response and cellular activity. In addition, we performed a correlation analysis of DEMs and differentially expressed genes in sea cucumber and 101 key negative regulator genes of DEMs were obtained, which may play important roles in regulating the process of A. japonicus under heat stress. Our study provided an increasing understanding of miRNAs’ roles in A. japonicus during heat stress. Through the discovery of their related targets, we have a deeper understanding of miRNAs and functional genes. In conclusion, our results provided new insights into the miRNA regulation and molecular adaptation mechanisms of A. japonicus under heat stress.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Sequence read archive under the accession number SRR18918495 ∼ SRR18918500.

Author contributions

MC, BL, ML, and XR contributed to conception and design of the research. MC and BL performed experiment. MC, BL, YW, XR, and JW performed data processing and statistical analysis. MC, ZZ, CW, and YY drafted the manuscript. XR and ML contributed to revision of manuscript for important intellectual content. ML and XR gave laboratory and project support. And all authors read and approved the final manuscript.

Funding

This work was funded by the National Key R&D Program of China (2018YFD0900305), Agriculture Seed Improvement Project of Shandong Province (2020LZGC015) and Key R&D Program of Qingdao (22-3-3-hygg-1-hy) and Central Public-Interest Scientific Institution Basal Research Fund, CAFS (2020TD40 and 2021GH05).

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.

Supplementary material

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

Supplementary Figure S1 | PCA analysis of six body wall libraries in Apostichopus japonicas.

Supplementary Figure S2 | Length distribution of valid reads in six libraries in the body wall of Apostichopus japonicas.

Supplementary Figure S3 | Families and members distribution of known miRNAs.

References

Bartel D. P. (2004). MicroRNAs. Cell. 116, 281–297. doi:10.1016/S0092-8674(04)00045-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Berthelot C., Brunet F., Chalopin D., Juanchich A., Bernard M., Noël B., et al. (2014). The rainbow trout genome provides novel insights into evolution after whole-genome duplication in vertebrates. Nat. Commun. 5, 3657. doi:10.1038/ncomms4657

PubMed Abstract | CrossRef Full Text | Google Scholar

Boonchuen P., Maralit B. A., Jaree P., Tassanakajon A., Somboonwiwat K. (2020). MicroRNA and mRNA interactions coordinate the immune response in non-lethal heat stressed Litopenaeus vannamei against AHPND-causing vibrio parahaemolyticus. Sci. Rep. 10, 787. doi:10.1038/s41598-019-57409-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan S. Y., Loscalzo J. (2010). MicroRNA-210: A unique and pleiotropic hypoxamir. Cell cycle. 9, 1072–1083. doi:10.4161/cc.9.6.11006

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen K., Cheng H.-H., Zhou R.-J. (2012). Molecular mechanisms and functions of autophagy and the ubiq-uitin-proteasome pathway. Hered. (Beijing) 34, 5–18. doi:10.3724/SP.J.1005.2012.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedländer M. R., Mackowiak S. D., Li N., Chen W., Rajewsky N. (2012). Mirdeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi:10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones-Rhoades M. W., Bartel D. P. (2004). Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol. Cell. 14, 787–799. doi:10.1016/j.molcel.2004.05.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang H., Liang Q.-J., Hu R., Li Z.-H., Liu Y., Wang W.-N. (2019). Integrative mRNA-miRNA interaction analysis associated with the immune response of Epinephelus coioddes to Vibrio alginolyticus infection. Fish Shellfish Immunol. 90, 404–412. doi:10.1016/j.fsi.2019.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Li C., Feng W., Qiu L., Xia C., Su X., Jin C., et al. (2012). Characterization of skin ulceration syndrome associated microRNAs in sea cucumber Apostichopus japonicus by deep sequencing. Fish shellfish Immunol. 33, 436–441. doi:10.1016/j.fsi.2012.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Li C., Zhao M., Zhang C., Zhang W., Zhao X., Duan X., et al. (2016b). miR210 modulates respiratory burst in Apostichopus japonicus coelomocytes via targeting Toll-like receptor. Dev. Comp. Immunol. 65, 377–381. doi:10.1016/j.dci.2016.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Li G., Ren L., Sun G., Yang J., Wei M. (2016a). Effects of hypoxic stress on oxidative stress indices in Apostichopus japonicus. Prog. Fish. Sci. 37, 133–139. doi:10.11758/yykxjz.20150708001

CrossRef Full Text | Google Scholar

Li X. (2021). Genomic DNA methylation level and transcriptome differences of sea cucumber under temperature and pathogen stress with their conjoint analysis. Shanghai: Shanghai Ocean University. [master’s thesis].

Google Scholar

Li X. (2020). “Study on pathogenicity of non-01 Vibrio cholerae to Macrobrachium nipponensis, host immune response and probiotic effect of antagonistic bacteria. Yangzhou: Yangzhou University. [doctor’s thesis].

Google Scholar

Li Y., Wang R., Li Y., Li Y., Mou C., Sun H., et al. (2018). Genome-wide profiling of DNA methylation in Apostichopus japonicus based on methylRAD-Seq. Periodical Ocean Univ. China Natural Sci. 48, 41–50. doi:10.16441/j.cnki.hdxb.20170242

CrossRef Full Text | Google Scholar

Livak K. J., Schmittgen T. D. (2013). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 25, 402–408. doi:10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma F., Liu Z., Huang J., Li Y., Kang Y., Liu X., et al. (2019). High-throughput sequencing reveals microRNAs in response to heat stress in the head kidney of rainbow trout (Oncorhynchus mykiss). Funct. Integr. Genomics. 19, 775–786. doi:10.1007/s10142-019-00682-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Manzerra P., Brown I. R. (1990). Time course of induction of a heat shock gene (hsp70) in the rabbit cerebellum after LSD in vivo: Involvement of drug-induced hyperthermia. Neurochem. Res. 15, 53–59. doi:10.1007/bf00969184

PubMed Abstract | CrossRef Full Text | Google Scholar

Ministry of Agriculture and Rural Affairs (2021). China fishery statistical yearbook. Beijing: China Agriculture Press.

Google Scholar

Naguibneva I., Ameyar-Zazoua M., Polesskaya A., Ait-Si-Ali S., Groisman R., Souidi M., et al. (2006). The microRNA miR-181 targets the homeobox protein hox-a11 during mammalian myoblast differentiation. Nat. Cell Biol. 8, 278–284. doi:10.1038/ncb1373

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohashi M., Holthaus A. M., Calderwood M. A., Lai C.-Y., Krastins B., Sarracino D., et al. (2015). The EBNA3 family of epstein-barr virus nuclear proteins associates with the USP46/USP12 deubiquitination complexes to regulate lymphoblastoid cell line growth. PLoS Pathog. 11, e1004822. doi:10.1371/journal.ppat.1004822

PubMed Abstract | CrossRef Full Text | Google Scholar

Olaru A. V., Selaru F. M., Mori Y., Vazquez C., David S., Paun B., et al. (2011). Dynamic changes in the expression of microRNA-31 during inflammatory bowel disease-associated neoplastic transformation. Inflamm. Bowel Dis. 17, 221–231. doi:10.1002/ibd.21359

PubMed Abstract | CrossRef Full Text | Google Scholar

Ou J., Li Y., Ding Z., Xiu Y., Wu T., Du J., et al. (2013). Transcriptome-wide identification and characterization of the Procambarus clarkii microRNAs potentially related to immunity against spiroplasma eriocheiris infection. Fish Shellfish Immunol. 35, 607–617. doi:10.1016/j.fsi.2013.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen I. M., Cheng G., Wieland S., Volinia S., Croce C. M., Chisari F. V., et al. (2007). Interferon modulation of cellular microRNAs as an antiviral mechanism. Nature. 449, 919–922. doi:10.1038/nature06205

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi J., Qiao Y., Wang P., Li S., Zhao W., Gao C. (2012). microRNA-210 negatively regulates LPS-induced production of proinflammatory cytokines by targeting NF-κB1 in murine macrophages. FEBS Lett. 586, 1201–1207. doi:10.1016/j.febslet.2012.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Reuss D. E., Piro R. M., Jones D. T. W., Simon M., Ketter R., Kool M., et al. (2013). Secretory meningiomas are defined by combined klf4 k409q and traf7 mutations. Acta Neuropathol. 125 (3), 351–358. doi:10.1007/s00401-013-1093-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosati A., Graziano V., De Laurenzi V., Pascale M., Turco M. C. (2011). BAG3: A multifaceted protein that regulates major cell pathways. Cell Death Dis. 2, e141. doi:10.1038/cddis.2011.24

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon P., Markiel A., Ozier O., Baliga N., Wang J., Ramage D., et al. (2003). Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome. Res. 11, 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao Y., Wang Z., Lv Z., Li C., Zhang W., Li Y., et al. (2016). NF-κB/Rel, Not STAT5, Regulates Nitric Oxide Synthase Transcription in Apostichopus japonicus. Dev. Comp. Immunol. 61, 42–47. doi:10.1016/j.dci.2016.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi T., Meng P., Qiang F., Shi H., Zhang H., Ren Y., et al. (20142014). Construction and preliminary identification of lentiviruses overexpressing and inhibiting Bos taurus miR-193a. Biotechnology. 24 52–57.

Google Scholar

Song H., Qi L., Zhang T., Wang H.-y. (2017). Understanding microRNA regulation involved in the metamorphosis of the veined rapa whelk (Rapana venosa). G3 Genes| Genomes |Genetics. 7, 3999–4008. doi:10.1534/g3.117.300210

CrossRef Full Text | Google Scholar

Suárez Y., Wang C., Manes T. D., Pober J. S. (2010). Cutting edge: TNF-induced microRNAs regulate TNF-induced expression of e-selectin and intercellular adhesion molecule-1 on human endothelial cells: Feedback control of inflammation. J. I. 184, 21–25. doi:10.4049/jimmunol.0902369

CrossRef Full Text | Google Scholar

Tian W., Wang G., Liu Y., Huang Z., Zhang C., Ning K., et al. (2017). Retracted: The miR-599 promotes non-small cell lung cancer cell invasion via SATB2. Biochem. Biophysical Res. Commun. 485, 35–40. doi:10.1016/j.bbrc.2017.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang T., Yang H., Zhao H., Chen M., Wang B. (2011). Transcriptional changes in epigenetic modifiers associated with gene silencing in the intestine of the sea cucumber, Apostichopus japonicus (Selenka), during aestivation. Chin. J. Ocean. Limnol. 29, 1267–1274. doi:10.1007/s00343-011-0143-2

CrossRef Full Text | Google Scholar

Wang X., Jin Y., Zhang H., Huang X., Zhang Y., Zhu J. (2018). MicroRNA-599 inhibits metastasis and epithelial-mesenchymal transition via targeting EIF5A2 in gastric cancer. Biomed. Pharmacother. 97, 473–480. doi:10.1016/j.biopha.2017.10.069

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang Y., Sui Y., Zhu Q., Sui X. (2017). Hsa-miR-599 suppresses the migration and invasion by targeting brd4 in breast cancer. Oncol. Lett. 14, 3455–3462. doi:10.3892/ol.2017.6651

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei M., Xu W.-t., Li K.-m., Chen Y.-d., Wang L., Meng L., et al. (2018). Cloning, characterization and functional analysis of dctn5 in immune response of Chinese tongue sole (Cynoglossus semilaevis). Fish Shellfish Immunol. 77, 392–401. doi:10.1016/j.fsi.2018.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen Z., Zou S., Chen M., Zhou H., Sun G., Feng Y., et al. (2021). DNA methylation level of genomic DNA of Apostichopus japonicus at different temperatures. Prog. Fish. Sci. 42, 46–54. doi:10.19663/j.issn2095-9869.20201204002

CrossRef Full Text | Google Scholar

Wienholds E., Kloosterman W. P., Miska E., Alvarez-Saavedra E., Berezikov E., de Bruijn E., et al. (2005). MicroRNA expression in zebrafish embryonic development. Science. 309, 310–311. doi:10.1126/science.1114519

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu L.-G., Li L.-Y., Shu H.-B. (2004). TRAF7 potentiates MEKK3-induced AP1 and CHOP activation and induces apoptosis. J. Biol. Chem. 279, 17278–17282. doi:10.1074/jbc.C400063200

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu P., Guo H., Wang H., Xie Y., Lee S. C., Liu M., et al. (2019). Identification and profiling of microRNAs responsive to cadmium toxicity in hepatopancreas of the freshwater crab Sinopotamon henanense. Hereditas. 156, 34. doi:10.1186/s41065-019-0110-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu L., Yu S., Xiong Y., Dui Y., Wang D., Huang W. (2021). Effects of miR-210 on proliferation and Apoptosis of dental pulp stem cells based on PI3K/AKT Pathway. Prog. Mod. Biomed. 21, 1212–1216. doi:10.13241/j.cnki.pmb.2021.07.003

CrossRef Full Text | Google Scholar

Zhang L., Jia D., Wei M., Cao Y., Jie M., Neurosurgery D. (2018). The miR-3188 promotes proliferation and invasion of glioma cell by targeting WDR20. J. Shanxi Med. Univ. 49, 802–808. doi:10.13753/j.issn.1007-6611.2018.07.011

CrossRef Full Text | Google Scholar

Zhang Y., Yu S., Liao M., Dong Y. (2022). Evaluation and prediction of extreme high temperature on sea cucumber (Apostichopus japonicus) pond aquaculture in China. J. Fish. Sci. China 29, 408. doi:10.12264/JFSC2021-0376

CrossRef Full Text | Google Scholar

Zhao G., Hou J., Xu G., Xiang A., Kang Y., Yan Y., et al. (2017). Cellular microRNA miR-10a-5p inhibits replication of porcine reproductive and respiratory syndrome virus by targeting the host factor signal recognition particle 14. J. General Virology. 98, 624–632. doi:10.1099/jgv.0.000708

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao Y., Chen M., Storey K., Sun L., Yang H. (2015). DNA Methylation Levels Analysis in Four Tissues of Sea Cucumber Apostichopus Japonicus Based on Fluorescence-Labeled Methylation-Sensitive Amplified Polymorphism (F-MSAP) During Aestivation. Comp. Biochem. Phys. B. 181, 26–32. doi:10.1016/j.cbpb.2014.11.001

CrossRef Full Text | Google Scholar

Zhong L., Zhang F., Zhai Y., Cao Y., Zhang S., Chang Y. (2015). Identification and comparative analysis of complement C3-associated microRNAs in immune response of Apostichopus japonicus by high-throughput sequencing. Sci. Rep. 5, 17763. doi:10.1038/srep17763

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou Q. (2014). The regulation of Bovine herpesvirus 5 encoded miR-10 on its target gene Us3. Wuhan: Huazhong Agriculture University. [doctor’s thesis].

Google Scholar

Zhou X., Chang Y., Zhan Y., Wang X., Lin K. (2018). Integrative mRNA-miRNA interaction analysis associate with immune response of sea cucumber Apostichopus japonicus based on transcriptome database. Fish Shellfish Immunol. 72, 69–76. doi:10.1016/j.fsi.2017.10.031

PubMed Abstract | CrossRef Full Text | Google Scholar