key: cord-0734157-1roubx7a
authors: Huo, Da; Sun, Lina; Storey, Kenneth B.; Zhang, Libin; Liu, Shilin; Sun, Jingchun; Yang, Hongsheng
title: The regulation mechanism of lncRNAs and mRNAs in sea cucumbers under global climate changes: Defense against thermal and hypoxic stresses
date: 2020-03-20
journal: Sci Total Environ
DOI: 10.1016/j.scitotenv.2019.136045
sha: 308b85a5c3d4cad71564cb5f7612c515f393af9d
doc_id: 734157
cord_uid: 1roubx7a

The aquatic environment can be greatly impacted by thermal and hypoxic stresses, particularly caused by intensified global warming. Hence, there is an urgency to understand the response mechanisms of marine organisms to adverse environment. Although long non-coding RNAs (lncRNAs) are involved in many biological processes, their roles in stress responses still remain unclear. Here, differentially expressed (DE) lncRNAs and mRNAs were identified as responses to environmental stresses in the economically important sea cucumber, Apostichopus japonicus, and their potential roles were explored. Based on a total of 159, 355 and 495 significantly upregulated genes and 230, 518 and 647 significantly downregulated genes identified in the thermal, hypoxic and combination thermal + hypoxic stress treatments, respectively, we constructed DE-lncRNA-mRNA coexpression networks. Among the networks, eight shared pairs were identified from the three treatments, and based on the connectivity degree, MSTRG.27265, MSTRG.19729 and MSTRG.95524 were shown to be crucial lncRNAs. Among all the significantly changed lncRNAs identified by RT-qPCR and sequencing data, binding sites were found in four other lncRNAs (MSTRG.34610, MSTRG.10941, MSTRG.81281 and MSTRG.93731) with Aja-miR-2013-3p, a key miRNA that responds to hypoxia in sea cucumbers. The hypoxia-inducible factor (HIF-1α) was also shown as the possible targeted mRNA of Aja-miR-2013-3p. As indicated by a dual-luciferase reporter assay system, “HIF-1α gene/Aja-miR-2013-3p/MSTRG.34610” network and the “HIF-1α gene/Aja-miR-2013-3p/MSTRG.10941” network may play important roles in sea cucumbers under environmental stresses. Moreover, environmental stress altered the expression of multiple lncRNAs and mRNAs, thus affecting various biological processes in A. japonicus, including immunity, energy metabolism and the cell cycle. At the molecular level, more comprehensive responses were elicited by the combined thermal/hypoxic stress treatment than by individual stresses alone in sea cucumbers. This study lays the groundwork for future research on molecular mechanisms of echinoderm responses to thermal and hypoxic stress in the context of global climate changes.

Due to global climate change, the natural water temperature has significantly increased, resulting in decreased dissolved oxygen levels in the past few decades (Dong et al., 2020) . Thus, in aquaculture ponds, high temperature (HT) and low dissolved oxygen (LO) conditions in aquatic systems may co-occur more frequently. Extremely high temperature and low dissolved oxygen would cause thermal stress and hypoxic stress, respectively. These environmental stresses could lead to yield loss and reduce the germplasm quality in the aquaculture industry (Huo et al., 2018) . The sea cucumber Apostichopus japonicus is a deposit feeder that participates in carbon cycling in marine ecosystems (Huo et al., 2019b) . In China, it is an economically important food species raised extensively in the mariculture; for example, the area, yield and quantity of seedlings in the country were 2.2 × 10 5 hm 2 , 2.1 × 10 8 kg and 52.8 billion by 2017, respectively (Ministry of Agriculture, 2018). Thermal and hypoxic stresses have been identified as the main causes of massive death of sea cucumbers over the summer in China in recent years, with these stresses negatively impacting of these stresses on the behavioral, physiological, histomorphological and molecular characteristics of A. japonicus (Huo et al., 2019a) . To reduce the damage caused by environmental stress, sophisticated adaptive response mechanisms have been developed by sea cucumbers. The regulation of gene expression occurs at the transcriptional, posttranscriptional and posttranslational levels (Shukla et al., 2008) .

Long non-coding RNAs (LncRNAs) are RNA molecules without protein-coding abilities and lengths over 200 nucleotides (Kowalczyk et al., 2012) . LncRNAs can function at the site of transcription in cis (cis-regulation) or trans (trans-regulation) modes and, as such, participate in multiple networks regulating gene expression through diverse mechanisms including alternative splicing (Gonzalez et al., 2015) , translation (Carrieri et al., 2012; Yoon et al., 2012) , and regulation of epigenetic and posttranscriptional gene (Ibeagha-Awemu et al., 2018; Mercer and Mattick, 2013; Yoon et al., 2013) . Although lncRNAs have been reported to play roles in the regulation of the innate immunity (Heward and Lindsay, 2014) , chromatin reprogramming (Gupta et al., 2010) , cellular differentiation and development (Fatica and Bozzoni, 2013) , as well as cardiac development and aging (The Cardiolinc et al., 2015) in model organisms, the biological functions of lncRNAs remain largely unknown in aquatic animals, especially for those involved in stress-responsive mechanisms.

Understanding dysregulated lncRNAs and mRNAs is important for revealing the molecular mechanisms occurring in sea cucumber coping with stress. Due to the recent availability of high-quality genomic information for A. japonicus , research on non-coding RNAs in sea cucumber can now be conducted. Therefore, we obtained lncRNA and messenger RNA (mRNA) expression profiles from respiratory tree of sea cucumbers under thermal stress, hypoxic stress, and combined thermal + hypoxic stress (HL) as well as control treatment conditions (NC) . Respiratory tree was selected as the target tissue in the present study, because it is the major organ responsible for respiratory metabolism in the sea cucumber and is very sensitive to environmental stress . The integrated analysis of the differentially expressed (DE) lncRNAs and mRNAs may provide clues to find genes with active roles in stress adaption in sea cucumbers under environmental stress in the context of global climate change. Our results represent the first construction of lncRNA-mRNA coexpression networks and lncRNA-miRNA-mRNA networks identified by the dualluciferase reporter system in sea cucumbers. The data reported herein provides important information about stress resistance in echinoderms with potential applications to other species in aquaculture.

In the present study, before using the sea cucumbers collected from coast of Weihai, China, one-week acclimation was proceeded in tanks containing aerated sand-filtered seawater at 16 ± 0.5°C. Around 50 active sea cucumbers (wet weights: 90-110 g) were then randomly selected and divided into groups. Sea cucumbers from each group were cultivated in separated tanks and exposed to stress for 48 h, and the moment when the water temperature had risen and/or oxygen decreased to the corresponding set value was defined as the treatment start time. For A. japonicus, water with sufficient dissolved oxygen (8 mg/L) and a temperature of 16°C were suitable for survival, thus this condition was taken as control (NC). Treatment group with low dissolved oxygen (2 mg/L) condition was set as the hypoxia treatment group (LO). The thermal stress treatment group (HT) was maintained at a high temperature of 26°C and the thermal and hypoxia combined treatment group (HL) was maintained at a high temperature of 26°C and low dissolved oxygen concentration (2 mg/L). Acclimation was carried out by using a 2-kW heating rod and dissolved oxygen control system, respectively (Huo et al., 2018) . Respiratory trees from nine sea cucumber individuals in each group were sampled promptly and preserved in liquid nitrogen, and stored at −80°C. The sea cucumber (A. japonicus) is not an endangered or a protected species and no permission was needed for sea cucumber collection.

Total RNA from tissues was extracted using MiniBEST Universal RNA Extraction Kit (Takara, Shiga, Japan) following the manufacturer's recommendations. RNA degradation and purity were detected using 1% agarose gels and kaiaoK5500® spectrophotometer (Kaiao, Beijing, China), respectively. The RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) was used for RNA integrity and concentration assessment. All samples had an RNA integrity number (RIN) ≥7.0, OD260/OD280 = 2.0 ± 0.2 and OD260/OD230 = 2.0 ± 0.2.

For the RNA sample preparations, 3 μg of RNA per sample was used. The sequencing libraries were constructed with varied index labels by a NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, USA) , and the details are as follows: first, ribosomal RNA was removed utilizing Epicentre Ribo-Zero™ Gold Kits (Epicentre, USA), and then NEBNext First-Strand Synthesis Reaction Buffer under raised temperature was used to fracture RNA into short RNA strands. Subsequently, random hexamer primers were used to synthesize the firststrand cDNA by utilizing RNA fragments as a template. Afterwards, second-strand cDNA synthesis was performed using a buffer, dNTPs, DNA polymerase I and RNase H. Then the library fragments were purified with QIAQuick PCR kits, and eluted with elution buffer (EB). Later on, terminal repair, poly(A) tailing and adapter ligation were implemented. The library fragments were purified with agarose gel electrophoresis to referentially select cDNA fragments of 300 bp in length and the UNG enzyme was used to digest second strand cDNA. After performing PCR, aimed products were retrieved by agarose gel electrophoresis. After measuring the RNA concentration of each library by a Qubit® RNA Assay Kit with a Qubit® 2.0 instrument, the RNA was diluted to 1 ng/μl. The insert size was assessed using an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA), and the Taqman fluorescence probe of the AB Step One Plus Real-Time PCR system (library valid concentration N 10 nM) was used to quantify the qualified insert size. Clustering of the index-coded samples was performed on a cBot cluster generation system using a TruSeq PE Cluster Kit v4-cBot-HS (Illumina). After generating cluster, the libraries were sequenced on an Illumina HiSeq 4000 platform and 150 bp paired-end reads were generated.

Raw data were processed with Perl scripts to ensure the quality of data to be used in subsequent analyses. Adaptor-polluted reads (reads containing N5 adapter-polluted bases), low-quality reads (reads with low-quality bases (phred quality value b19) accounting for N15% of the total bases) and the reads with N bases accounting for N5% of total bases were removed. For the paired-end sequencing data, both reads were filtered out if any read of the paired-end reads was adaptorpolluted. The filtered clean data was subjected to statistical analyses for quantity and quality.

The reference genomes and the annotation file of A. japonicus were used . Bowtie2 (v2.2.3) was used for building the genome index, and the clean data were mapped to the reference genome using TopHat (v2.0.12). Moreover, TopHat could call Bowtie2 for mapping, increasing speed and accuracy. Cufflinks was used to predict the new transcripts and those with lengths were N200 bp, exons N2 and cover degree ≥3 were selected. Coding Potential Calculator (CPC), Coding-Non-Coding-Index (CNCI), Coding Potential Assessment Tool (CPAT) and Pfam were used to identify lncRNAs, and those lncRNAs that were identified in only one sample were removed. Venn diagrams were presented using Venny 2.1 (Oliveros, 2007) .

Fragments per kilobase million (FPKM) were then calculated to estimate the expression level of genes in each sample. The formula is as follows:

where N is the total number of mapped fragments in the given sample, L is the length of the certain gene, and F is the number of fragments in a given sample that is assigned to a certain gene. FPKM can eliminate the effect of sequencing depth and gene length on the gene expression level, and the data can be compared between each other directly.

DESeq (v1.16) was used for differential gene expression analysis between two samples with biological replicates using a model based on the negative binomial distribution. A P-value was assigned to each gene and adjusted by Benjamini and Hochberg's approach for controlling the false discovery rate as a Q-value. Genes with q b 0.05 and | log2FoldChange | ≥ 1 were identified as differentially expressed genes.

The GO (http://www.geneontology.org/) and KEGG (http://www. genome.jp/kegg/) databases were used to classify and group the identified mRNAs and lncRNAs. GO and KEGG pathway enrichment analyses were based on significantly up-and down-regulated genes, respectively. The ratio of the number of differentially regulated genes annotated in a pathway term to the number of all genes annotated in the same pathway term was defined as rich factor. A greater rich factor indicates greater intensiveness. The hypergeometric test was used to identify the significantly enriched GO terms and pathways of differentially regulated genes, and those with P-value b 0.05 were considered significant (Huo et al., 2019b) . The GO terms with top 5 lowest Pvalues for the ontologies biological process, molecular function and cellular component in the three treatments showing in Fig. 3 were made by ggplot2 (Wickham, 2016) . The most significantly enriched GO terms and most significant KEGG pathway were identified with the lowest P-value among all the significantly enriched terms.

Ten DE-mRNAs and ten DE-lncRNAs with high-fold-change values in the three group comparisons were selected for expression detection. Total RNA was extracted from the respiratory tree, the same tissue used for the construction of the RNA-seq profile, by using a MiniBEST Universal RNA Extraction Kit (Takara, Shiga, Japan). A NanoDrop 1000 (Thermo Fisher Scientific, USA) was used for measuring the quality and concentration of RNA before using reverse transcriptase (Takara) to synthesize first-strand cDNA. A SYBR Green® real-time PCR assay (SYBR PrimeScript ™ RT-PCR Kit II, TaKaRa) with an Eppendorf Mastercycler® ep realplex (Eppendorf, Hamburg, Germany) were used for examining the mRNA expression levels, and β-actin was used as a reference gene for internal standardization (Yang et al., 2010b) .

Primer3 was used for designing primers according to the sequence information in the transcriptome database. All primer sequences are shown in Table S1 . In the present study, three biological replicates × three technical replicates were used for the detection of gene expression in each group. The total volume of the amplification mix was 20 μL, which contained 10 μL of SYBR Green Master Mix (Takara), 8 μL of RNase-free water, 0.5 μL (each) of the forward and reverse primers (10 mM), and 1 μL of diluted cDNA. The thermal cycling was performed according to the following procedure: 95°C for 5 s, followed by 40 cycles at 95°C for 10 s, 60°C for 20 s and 72°C for 30 s. Melting curve analysis was used for illustrating the specificity of the amplification products, and the 2 −△△CT method was used to analyze the comparative mRNA expression levels of the selected genes (Schmittgen and Livak, 2008) . All data were analyzed by SPSS19 software (IBM Corp., Armonk, NY, USA), and presented as mean ± SD. The threshold of statistical significance was a P-value b 0.05 analyzed by one-way ANOVA with Tukey's test. Among the RT-qPCR results of the three treatments compared with the control, genes showing the same tendency as the sequencing data for all the three stress treatments relative to the control treatment were defined as completely matched, and those showing the same tendency in at least one comparison were defined as partially matched.

To identify the trans-targets of lncRNAs, the FASTA sequences of lncRNAs and mRNAs were extracted based on the genome and annotation file, and they were used to calculate the free energy of binding using RNAhybrid, a software specially designed to quickly find possible hybridization sites for a query RNA in large RNA databases and to predict the RNA-RNA interaction (Rehmsmeier et al., 2004) . The lncRNA-mRNA co-expression network was built by Cytoscape software (Cytoscape Consortium, San Diego, CA, USA) (Shannon et al., 2003) based on the correlation between the DE-lncRNAs and DE-mRNAs. LncRNAs and coding genes were identified based on Pearson's correlation coefficients (equal to or N0.8) for trans-targets or based on the genomic positional relation in 50 kb regions (upstream and downstream) for cis-targets. Based on our previous study, Aja-miR-2013-3p was taken as a key stress-responsive-miRNA connected with lncRNAs and mRNAs in A. japonicus . Hypoxia-inducible factor (HIF-1α), a key hypoxia responder, was predicted by RNAhybrid to be one of the targets of Aja-miR-2013-3p. Thus, these two were used for lncRNA binding possibility screening and lncRNA-miRNA-mRNA network prediction. All the plasmids used in the present study were custom synthesized by GeneChem Co., Ltd. (Shanghai, China) (Table S2) and cloned into a GV272 firefly luciferase plasmid (GeneChem). The GV272 empty vector was utilized as a negative control (3′ UTR-NC) (GeneChem), and the multiple cloning sites of it were digested by XbaI/XbaI restriction enzymes. The 3′ untranslated regions (UTR) segment of mRNA (HIF-1α) and lncRNA (MSTRG.93731, MSTRG.10941, MSTRG.34610 and MSTRG.81281) or their mutant was inserted into the vector to yield the 3′ UTR or 3′ UTR-MU. The recombinant products were verified by DNA sequencing, and the group details can be found in Table S2 . Hsa-miR-146b and its target gene TRAF6 were utilized as positive controls (Park et al., 2015) . The pGL3 control vector (Promega) and pRL-TK were used as a firefly luciferase reporter vector and a Renilla luciferase control reporter vector, respectively. PCR amplification of the RNA sequences was performed using a PrimerSTAR Max DNA Polymerase Mix (Takara Bio Inc., Japan). The top 10 competent bacterial cells were transformed with the ligation product of the luciferase reporter vector and the digested PCR fragment. Transformed cells were incubated at 37°C on LB agar plates with antibiotics overnight, and isolated clones were inoculated in LB medium with antibiotics at 37°C with shaking overnight. The ZR Plasmid Miniprep Kit (Zymo Research, Orange, CA, USA) was used for the plasmid DNA purify. All constructs were verified by sequencing (Table S2 ). The transfection and detection of the luciferase activities were performed using Li's method . The human embryonic kidney cell line, HEK 293 T at 90% confluence in 24-well plates was transfected with a 3′ UTR luciferase plasmid (0.1 μg), a microRNA plasmid (0.4 μg) and a Renilla plasmid (0.02 μg) using the X-tremeGENE HP DNA Transfection Reagent (Roche, Basel, Switzerland). The firefly luciferase plasmid and the Aja-miR-2013-3p expression plasmid were co-transfected with a Renilla luciferase vector (Promega, Madison, WI, USA) for normalization. Luciferase activity was measured using a microplate reader (M2009PR, Tecan infinite) after 48 h following the manufacturer's instructions. All experiments were performed in triplicate, and the original data and the normalization method are shown in Table S2 . The analysis methods were the same as those desribed for RT-qPCR, including software, statistical test and figure construction.

3.1. RNA-sequencing and identification of expressed mRNA and lncRNA genes in the respiratory tree of A. japonicus

In the present study, RNA-seq of A. japonicus samples from control (NC), high temperature exposure (HT), low oxygen exposure (LO), and high temperature + low oxygen (HL) was performed, and the data were submitted to the GEO database [#GSE131676]; 118 to 159 million raw reads and 105 to 143 million clean reads were obtained per sample. Moreover, 28,586 mRNAs were identified. Every potential lncRNA in all samples was assembled by Stringtie. Finally, 59,967 lncRNAs were retained at the intersection of the CPC, CNCI, CPAT and Pfam scans (Fig. S1 ). Because of the lack of lncRNA information for invertebrates, especially for sea cucumbers, all identified lncRNAs were novel, and the names of novel lncRNAs begin with "MSTRG". Table S3 . Key features of DE-lncRNAs, such as exon number and sequence length, were characterized by comparison with identified DE-mRNAs under the same conditions (Fig. 1) . The results showed that the exon number of lncRNAs was less than that of mRNAs. Moreover, most lncRNAs had fewer than 7 exons, but several mRNAs had more than thirty exons. In addition, the exon numbers of lncRNAs were distributed over a wider range than mRNAs. The lengths of most lncRNAs were shorter than mRNAs.

A total of 389, 873 and 1142 representative DEGs (mRNAs and lncRNAs) were identified in sea cucumbers under thermal, hypoxic and thermal + hypoxic stress compared with controls based on the criteria of |log 2 FoldChange| ≥ 1 and q b 0.05. When A. japonicus was exposed to thermal stress, 88 mRNAs and 71 lncRNAs were significantly upregulated and 230 genes were significantly downregulated, including 107 mRNAs and 123 lncRNAs. When A. japonicus was exposed to hypoxic stress, 355 genes were significantly upregulated including 247 mRNAs and 108 lncRNAs, whereas 518 were significantly downregulated, including 233 mRNAs and 285 lncRNAs. When A. japonicus was exposed to both thermal and hypoxic stress, more significantly changed genes were observed: 495 significantly upregulated genes were identified including 277 mRNAs and 218 lncRNAs, as well as 647 significantly downregulated genes including 346 mRNAs and 301 lncRNAs. Volcano plots illustrating these data are shown in Fig. S3 . The top 10 DE-lncRNAs and DE-mRNAs with highest fold change value are shown in Tables 1 and 2.

To confirm the expression of lncRNAs and mRNAs, RT-qPCR analysis was applied to verify the RNAseq results (Fig. 2) . Among the 10 selected lncRNAs and 10 mRNAs (Tables 1 and 2) , the expression levels of six lncRNAs and seven mRNAs completely matched those of the high-throughput sequencing data; these were MSTRG.10941, MSTRG.2005, MSTRG.83990, MSTRG.34610, MSTRG.93731, MSTRG.72169, complement factor B (AJAP22103), egg coat matrix protein (AJAP09945), serine/threonine-protein kinase Nek4 isoform X5 (AJAP24609), G2/mitotic-specific cyclin-B-like (AJAP03163), fatty acid-binding protein 2 (AJAP18992), topoisomerase IIassociated protein PAT1 (AJAP03092) and twin BRCT domain (AJAP10472). The expression levels of three other lncRNAs and one mRNA partially matched those of high-throughput sequencing data (MSTRG.79659, MSTRG.39061, MSTRG.26009, and heat shock protein 70). The expression levels of one lncRNA and two mRNAs did not achieve significant levels. Validation of key DE-lncRNA and DE-mRNA transcripts by RT-qPCR yielded results consistent with the RNA-seq analysis in most cases, thus confirming the reliability of the RNA-seq technique and the high quality of the identified lncRNAs and mRNAs.

GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted to explore the function of the DE-lncRNAs and DE-mRNAs ( Fig. 3 and Table 3 ). From the perspective of DE-mRNAs in the HT, LO and HL groups, the most significantly enriched GO terms in biological processes category were the cell cycle process, glucosaminecontaining compound biosynthetic process and the mitotic cell cycle; the most significantly enriched GO terms in cellular components category based on DE-mRNAs were sarcoplasmic reticulum lumen, extracellular region part and extracellular region; finally, the most significantly enriched GO terms in molecular functions category based on DE-mRNAs were metalloendopeptidase activity, monosaccharide binding and peptidase regulator activity.

Based on the results of significant DE-lncRNA analysis in the HT, LO and HL groups, the most significantly enriched GO terms in biological processes category were the regulation of Fc receptor mediated stimulatory signaling pathway, ribosomal large subunit biogenesis and microtubule-based process; the most significantly enriched GO terms in cellular components category based on DE-lncRNAs were ribosome, preribosome and intracellular non-membrane-bounded organelle. The most significantly enriched GO terms in molecular functions category based on DE-lncRNAs in the HT and HL groups was motor activity, and whereas that in the LO group was AMP binding.

Regarding KEGG pathways, based on the DE-mRNAs, olfactory transduction was the most significant pathway in the HT group, and thyroid hormone synthesis was the most significant pathway in LO and HL groups. Based on the DE-lncRNAs, ribosome was the most significant pathway in the HT and HL groups, and ribosome biogenesis in eukaryotes was the most significant pathway in the LO group (Table 3) .

3.4.1. LncRNA target gene prediction and lncRNA-mRNA co-expression network construction

In the present study, the cis-and trans-targets of the DE-lncRNAs were predicted to address the function of lncRNAs in concert with their target genes transcripts (mRNAs) to adapt to adverse environments (Table S4) . Among the trans-target gene predictions in the three groups, the number of positive correlations between the lncRNA-mRNA pairs was higher than the number of negative correlations. DE-lncRNAs and their corresponding DE cis-and trans-target genes were used to construct a lncRNA-mRNA co-expression network (Fig. 4) . The network included 42 lncRNA probes and 58 mRNA probes in the HT group, 178 lncRNA probes and 278 mRNA probes in the LO group, and 252 lncRNA probes and 405 mRNA probes in the HL group. The number of DE pairs found in the HT, LO and HL groups was 87, 1326 and 2430, respectively (Table S5 ). The greatest number of DE pairs in the HL group might be a result of the highest number of DEGs identified in the HL group. These results also indicated that the combined stress (HL) caused more extensive responses at the molecular level in sea cucumbers than the single stresses alone (HT and LO). Moreover, eight co-identified DE lncRNA-mRNA pairs were found in all three treatments (Fig. 4d) . Based on the connectivity degree, the most important lncRNAs in the interaction networks in the HT, LO and HL groups were 

In the present study, a dual-luciferase reporter assay system (DLRAS) was used for validation, and the results are shown in Fig. 5 . The group interactions without practical meanings or significance were not shown in the figure. Our previous study found several environmental-stress-responsive miRNAs that include Aja-miRNA- 2013-3p , and its predicted binding structure and mature sequence were shown in Fig. 6c . We analyzed the possible targeted mRNAs and lncRNAs, and found that hypoxia-inducible factor 1α (HIF-1α) was a potential target gene, which is a key responder in hypoxia defense. Moreover, binding sites for Aja-miRNA-2013-3p were also found in four lncRNAs (MSTRG.93731, MSTRG.10941, MSTRG.34610 and MSTRG.81281) (Fig. 6b) . Details of information of DLRAS, such as group content, sequences and primary data can be found in Table S2 . Compared with that in the positive control (PC)-miRNA-NC group, the PC-miRNA was significantly downregulated. Thus, the transfection system used in the current study proved to be effective. Compared with the mmi2 group, the mmi4 group showed a significant decrease of over 20% in expression. Thus, Aja-miRNA-2013-3p was proven to suppress the expression of HIF-1α. Compared with the lmi2 group, the lmi4 and lmi16 groups showed significant decreases of N20% in expression. Moreover, the lmi12 group showed significantly decreased expression levels compared with the expression levels in the lmi2 group, but the difference was only approximately 17%. Thus, Aja-miRNA-2013-3p was proven to bind lncRNAs such as MSTRG.93731, MSTRG.10941 and MSTRG.34610. Under hypoxic stress, the HIF-1α gene and MSTRG.10941 showed increased expression levels and could both bind Aja-miR-2013-3p; under thermal stress and combined stress, the HIF-1α gene and MSTRG.34610 showed decreased expression levels and could both bind Aja-miR-2013-3p. Therefore, the "HIF-1α gene/ Aja-miR-2013-3p/MSTRG.34610" network and the "HIF-1α gene/Aja-miR-2013-3p/MSTRG.10941" network might play important roles in sea cucumbers under environmental stress (Fig. 6) . The relative expression levels of these lncRNAs, miRNAs and mRNAs were compared for the three treatment groups and the NC group, and the results are shown in Fig. S4 (unpublished data).

Long noncoding RNAs (lncRNAs) could act as competing endogenous RNAs (ceRNAs) that compete with mRNAs for binding to miRNAs by sharing at least one miRNA response element (MRE), and thus affecting gene expression (Salmena et al., 2011) . Emerging evidence indicates that lncRNAs play crucial roles in many biological processes and contribute to stress resistance. However, the regulatory mechanisms remain unclear, especially in sea cucumber. To date, many lncRNAs have been identified in model animals including human (Wapinski and Chang, 2011) , rat (Wang et al., 2014) , zebrafish (Pauli et al., 2012) , arabidopsis (Liu et al., 2012) , fruit fly (Nyberg and Machado, 2016) , and elegans (Nam and Bartel, 2012) ; they have also been identified in other vertebrates, such as pig (Wang et al., 2016) and chicken (Li et al., 2012) ; and other marine animals, such as oyster (Feng et al., 2018) and salmon (Valenzuela-Miranda and Gallardo-Escárate, 2016). However, information concerning the systematic identification and function of lncRNAs in marine invertebrates remains lacking, especially for echinoderms. A. japonicus is a good model holothurian species, because of its economic and ecological value. LncRNAs of A. japonicus have fewer exons and shorter transcripts in length compared with mRNAs, and these results are in agreement with those of previous studies examining pig (Wang et al., 2016) , rat (Wang et al., 2014) and oyster (Feng et al., 2018) . With the absence of lncRNA data for other holothurians, the conservation of lncRNAs is difficult to properly evaluate in sea cucumbers. Delightfully, recent advances in sea cucumber genome sequencing have resulted in the identification of novel lncRNAs in A. japonicus and the exploration of novel regulatory pathways in different situations .

In the current study, based on sequencing data, we investigated transcriptomic changes and systematically identified the lncRNAs associated with the responses to environmental stress, specially, thermal, hypoxic and combined stresses. Although the functions of most lncRNAs in A. japonicus are not yet clear, evidence suggests that lncRNAs can regulate protein-coding gene expression in organisms through cis-acting mechanisms or trans-acting mechanisms (Han et al., 2012) . Thus, we constructed a coexpression network joint by lncRNAs and their target mRNAs and explored potential lncRNA-miRNA-mRNA networks, which may play crucial roles in A. japonicus exposed to adverse environmental stress. As a key modulator of most oxygen-limitation responses, HIF-1α has received extensive attention in hypoxiarelated studies. Aja-miRNA-2013-3p was reported to be important in sea cucumbers coping with hypoxic stress . The luciferase assay results showed that Aja-miRNA-2013-3p was able to bind to the HIF-1α gene, and so were lncRNAs, such as MSTRG.93731, MSTRG.10941 and MSTRG.34610. For example, under hypoxic conditions, Aja-miRNA-2013-3p was downregulated, but the HIF-1α gene and MSTRG.10941 were upregulated under hypoxic stress . This might be explained by a competent relationship between the HIF-1α gene and MSTRG.10941 binding to Aja-miRNA-2013-3p. However, no reported function was found for the novel lncRNAs mentioned above. Future studies of the regulatory mechanisms of coding and non-coding genes are needed.

Thermal stress and hypoxic stress have been shown to impact the immune response of organisms (Baze et al., 2011; Jin et al., 2011) . In vertebrates, emerging evidence suggests that lncRNAs are expressed in a highly lineage-specific manner and play an important role in Table 3 KEGG pathways analysis of differentially expressed mRNAs and lncRNAs in sea cucumbers under environmental stresses. immune regulation (Chen et al., 2017 ). The complement system is an essential immune response in invertebrates (Mastellos and Lambris, 2002) . In the present study, thermal stress and hypoxic stress lead to the up-regulation of complement factor B, indicating mediation of the alternative pathway. In addition, the complement component C3 was also found to be induced in sea cucumbers under environmental stress in our previous study (Huo et al., 2019b) . Heat shock proteins (HSPs) are widely viewed as a bioindicator of thermal stress in A. japonicus. HSPs can stimulate the immune system (Yang et al., 2010a) and were reported to be related to apoptosis (Roberts et al., 2010) and protein degradation (Nelson et al., 1992) . The expression of HSPs was more sensitive to thermal stress than to hypoxic stress, possibly because it is easier to induce protein denaturation during thermal stress than during hypoxic stress.

Thermal and hypoxic stress alters the body's energy metabolism because of trade-off effects. In the current study, many energy metabolism-related pathways, such as "amino sugar and nucleotide sugar metabolism", "fructose and mannose metabolism" and "fat digestion and absorption" were significantly changed under stress conditions (Table 3) . Fatty acid-binding protein 2 (FABP2) is a lipidbinding protein that mediates fatty acid absorption, transport and intracellular metabolism (Lowe et al., 1987; Nakanishi et al., 2004) .

Additionally, it is an important regulator of insulin resistance (Sipilainen et al., 1997) . Its significantly decreased expression under environmental stress (Fig. 2) also indicated the depression of fatty acid metabolism and may also indicate a potential changes of energy supply.

The cell cycle pathway was significantly impacted as indicated by KEGG analysis based on DE-mRNAs in the HT group and DE-lncRNAs in the HL group ( Fig. 3 and Table 3 ). Moreover, the cell cycle process was the most significantly enriched GO term of the biological processes category, as indicated by the DE-mRNAs in the HT, LO and HL groups. In the current study, topoisomerase II-associated protein (PAT1), a gene that accurately facilitates chromosome separation during cell division (Alzan et al., 2016) , showed significantly decreased expression under the three different environmental stresses relative its expression under control conditions (Fig. 2 ). In addition, thermal stress and the combined stress treatment significantly changed serine/threonine-protein kinase Nek4 isoform X5 (Nek4) expression, which potentially affects threonine residues, and is required for normal cell cycle arrest in response to double-stranded DNA damage (Nguyen et al., 2012) . All the evidence showed that the cell cycle process was extensively changed when sea cucumbers were exposed to environmental stress. 

Comparative bioinformatics analysis of transcription factor genes indicates conservation of key regulatory domains among Babesia bovis, Babesia microti, and Theileria equi

Chronic hypoxia stimulates an enhanced response to immune challenge without evidence of an energetic tradeoff

Long non-coding antisense RNA controls Uchl1 translation through an embedded SINEB2 repeat

Gene regulation in the immune system by long noncoding RNAs

Dietary supplementation of teprenone potentiates thermal and hypoxia tolerance as well as cellular stress protection of Epinephelus coioides juveniles reared under multiple stressors

Long non-coding RNAs: new players in cell differentiation and development

Transcriptional profiling of long non-coding RNAs in mantle of Crassostrea gigas and their association with shell pigmentation

A lncRNA regulates alternative splicing via establishment of a splicing-specific chromatin signature

Long non-coding RNA HOTAIR reprograms chromatin state to promote cancer metastasis

LncRNA profile of glioblastoma reveals the potential role of lncRNAs in contributing to glioblastoma pathogenesis

Long non-coding RNAs in the regulation of the immune response

Differential expression of miRNAs in the respiratory tree of the sea cucumber (Apostichopus japonicus) under hypoxia stress

Impact of hypoxia stress on the physiological responses of sea cucumber Apostichopus japonicus: respiration, digestion, immunity and oxidative damage

Metabolome responses of the sea cucumber Apostichopus japonicus to multiple environmental stresses: heat and hypoxia

Global-warming-caused changes of temperature and oxygen alter the proteomic profile of sea cucumber Apostichopus japonicus

Integration of lncRNA and mRNA Transcriptome analyses reveals genes and pathways potentially involved in calf intestinal growth and development during the early weeks of life

Chronic heat stress weakened the innate immunity and increased the virulence of highly pathogenic avian influenza virus H5N1 in mice

RNA discrimination

Identification of long non-protein coding RNAs in chicken skeletal muscle using next generation sequencing

Identification of circRNAs for miRNA targets by argonaute2 RNA immunoprecipitation and luciferase screening assays

Genome-wide analysis uncovers regulation of long intergenic noncoding RNAs in Arabidopsis

Expression of rat intestinal fatty acid-binding protein in Escherichia coli. Purification and comparison of ligand binding characteristics with that of Escherichia coli-derived rat liver fatty acidbinding protein

Complement: more than a 'guard' against invading pathogens?

Structure and function of long noncoding RNAs in epigenetic regulation

The effect of polymorphism in the intestinal fatty acid-binding protein 2 gene on fat metabolism is associated with gender and obesity amongst non-diabetic Japanese-Americans

Long noncoding RNAs in C. elegans

The translation machinery and 70 kd heat shock protein cooperate in protein synthesis

Nek4 regulates entry into replicative senescence and the response to DNA damage in human fibroblasts

Comparative expression dynamics of intergenic long noncoding RNAs in the genus Drosophila

VENNY. An Interactive Tool for Comparing Lists With Venn Diagrams

MicroRNA-146a and microRNA-146b regulate human dendritic cell apoptosis and cytokine production by targeting TRAF6 and IRAK1 proteins

Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis

Fast and effective prediction of microRNA/target duplexes

Heat shock proteins (chaperones) in fish and shellfish and their potential role in relation to fish health: a review

A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language?

Analyzing real-time PCR data by the comparative CT method

Cytoscape: a software environment for integrated models of biomolecular interaction networks

The role of microRNAs and other endogenous small RNAs in plant stress responses

Variants in the human intestinal fatty acid binding protein 2 gene in obese subjects

Long noncoding RNAs in cardiac development and ageing

Novel insights into the response of Atlantic salmon (Salmo salar) to Piscirickettsia salmonis: interplay of coding genes and lncRNAs during bacterial infection

Characteristics of long non-coding RNAs in the Brown Norway rat and alterations in the Dahl salt-sensitive rat

Analyses of long non-coding RNA and mRNA profiling using RNA sequencing during the preimplantation phases in pig endometrium

Long noncoding RNAs and human disease

ggplot2: Elegant Graphics for Data Analysis

Stability comparison of cytb and β-actin gene expression in sea cucumber Apostichopus japonicus

Expression of immune-related genes in embryos and larvae of sea cucumber Apostichopus japonicus

LincRNA-p21 suppresses target mRNA translation

Posttranscriptional gene regulation by long noncoding RNA

The sea cucumber genome provides insights into morphological evolution and visceral regeneration

The sea cucumber (A. japonicus) is not an endangered or a protected species and no permission was needed for sea cucumber collection.

The data generated during the current study were submitted to the GEO database [#GSE131676].

Although lncRNAs appear increasingly important in molecular mechanism studies, an understanding of lncRNAs in the echinoderm phylum remains lacking at present. To our knowledge, the present study represents the first comprehensive analysis of the lncRNA-mRNA co-expression networks and the validation of lncRNA-miRNA-mRNA networks in the response of marine organisms towards individual and combined influences of heat and hypoxia. Although this study was still limited by the absence of a lncRNA database for marine organisms, the identification and annotation of the putative lncRNAs provide a basis for the further clarification of the mechanisms underlying the regulation of mRNAs expression by lncRNAs in the environmental stress response of invertebrates. In future research, overexpression, CRISPR-Cas9 and RNA interference gene silencing strategies in sea cucumbers are expected to be useful for elucidating the specific roles of lncRNAs.Supplementary data to this article can be found online at https://doi. org/10.1016/j.scitotenv.2019.136045. 

Da Huo designed the study, carried out the experiment, analyzed data and drafted the paper.Lina Sun revised the draft and offered the foundation. Kenneth B. Storey contributed to examination of data and revision of the manuscript draft.Libin Zhang was responsible for the source of sea cucumber and offered the foundation.Shilin Liu contributed to the source of sea cucumber and sample preparation.Jingchun Sun helped with sequencing work and data curation. Hongsheng Yang designed the study and supervised the project.All authors contributed to manuscript revision, and read and approved the submitted version. 

The authors declare no conflict of interest.