- Research article
- Open Access
Transcriptomic analysis of sea cucumber (Holothuria leucospilota) coelomocytes revealed the echinoderm cytokine response during immune challenge
BMC Genomics volume 21, Article number: 306 (2020)
The sea cucumber Holothuria leucospilota belongs to echinoderm, which is evolutionally the most primitive group of deuterostomes. Sea cucumber has a cavity between its digestive tract and the body wall that is filled with fluid and suspended coelomic cells similar to blood cells. The humoral immune response of the sea cucumber is based on the secretion of various immune factors from coelomocytes into the coelomic cavity. The aim of this study is to lay out a foundation for the immune mechanisms in echinoderms and their origins in chordates by using RNA-seq.
Sea cucumber primary coelomocytes were isolated from healthy H. leucospilota and incubated with lipopolysaccharide (LPS, 10 μg/ml), polyinosinic-polycytidylic acid [Poly (I:C), 10 μg/ml] and heat-inactived Vibrio harveyi (107 cell/ml) for 24 h, respectively. After high-throughput mRNA sequencing on an Illumina HiSeq2500, a de novo transcriptome was assembled and the Unigenes were annotated. Thirteen differentially expressed genes (DEGs) were selected randomly from our data and subsequently verified by using RT-qPCR. The results of RT-qPCR were consistent with those of the RNA-seq (R2 = 0.61). The top 10 significantly enriched signaling pathways and immune-related pathways of the common and unique DEGs were screened from the transcriptome data. Twenty-one cytokine candidate DEGs were identified, which belong to 4 cytokine families, namely, BCL/CLL, EPRF1, IL-17 and TSP/TPO. Gene expression in response to LPS dose-increased treatment (0, 10, 20 and 50 μg/ml) showed that IL-17 family cytokines were significantly upregulated after 10 μg/ml LPS challenge for 24 h.
A de novo transcriptome was sequenced and assembled to generate the gene expression profiling across the sea cucumber coelomocytes treated with LPS, Poly (I:C) and V. harveyi. The cytokine genes identified in DEGs could be classified into 4 cytokine families, in which the expression of IL-17 family cytokines was most significantly induced after 10 μg/ml LPS challenge for 24 h. Our findings have laid the foundation not only for the research of molecular mechanisms related to the immune response in echinoderms but also for their origins in chordates, particularly in higher vertebrates.
Holothuria leucospilota is a tropical sea cucumber that belongs to phylum Echinodermata, class Holothuroidea, order Aspidochirotida and family Holothuriidae. Naturally, H. leucospilota is distributed in the Indo-West Pacific, mostly from eastern Africa to the Hawaii islands and Society islands in the Pacific ocean, and from southern Japan to the Sark bay in Australia . In recent years, H. leucospilota has become an emerging aquaculture species in southern China . Evolutionally, echinoderms are positioned taxonomically at the base of deuterostomes, along with the higher-order hemichordate and chordate groups . Therefore, studies on the biological processes in echinoderms, such as development, reproduction, metabolism and immunity, may provide new insights not only for echinoderms themselves, but also for the origins of these biological processes in chordates, particularly in higher vertebrates.
Given that lacking of the adaptive immunity, the innate immunity is the major mechanism for sea cucumber to defend the environmental pathogens. The innate immunity of sea cucumber includes multiple immune-related factors, such as, antimicrobial peptides (e.g. lectins, lysozyme, clotting protein and complement) [4, 5], antimicrobial reactive oxygen species , pattern recognition receptors, apoptosis [7, 8] and immune cytokines . Sea cucumber has a cavity between its digestive tract and body wall that is filled with coelomic fluid and suspended coelomocytes that are similar to the hematocyts of vertebrates. In sea cucumber, the cellular immunity is executed by coelomocytes, and the humoral immunity is based on a variety of macromolecules in the coelomic cavity that secreted by coelomocytes [10,11,12]. When sea cucumbers are infected by pathogenic microbes, they rely on their cellular and humoral immune responses to identify and eliminate the invading microbes and repair the wounds .
Next-generation sequencing (NGS) technology is a revolutionary change to traditional sequencing technology. The high-throughput mRNA sequencing (RNA-seq) is a transcriptomic research method that provides information on transcript expression and has the advantages of being easy-handle and low-cost . For non-model organisms, RNA-Seq is not limited to detecting transcripts that correspond to existing genomic sequence. RNA-Seq focus more on the coding region of genes and has very low background signal . Therefore, RNA-Seq is particularly attractive for non-model organisms with genomic sequences that are yet to be determined.
To date, transcriptomic sequencing technology has been widely applied for analyzing the gene expression profiles of echinoderms in a variety of developmental and physiological processes, and under multiple infected- or stressed-conditions. In the sea urchin Heliocidaris erythrogramma, RNA-seq has been used to measure the mRNA expression profiles of larvae, metamorphism and post-metamorphism life cycle stages, to elucidate the evolutional and developmental mechanism of the radial adult body in echinoderms . The expression dynamics across development has also been measured by RNA-seq in three sea urchin species: the lecithotroph H. erythrogramma, the closely related planktotroph H. tuberculata, and an outgroup planktotroph Lytechinus variegatus, to reveal how evolutionary changes in gene regulation contribute to phenotypic differences between different species . In addition, RNA-seq has been used to analyze the differentially expressed genes (DEGs) in the body cavity cells between the viral-infected and normal starfish Pycnopodia helianthoides, to illustrate the immune and nervous responses to the sea star wasting disease . Furthermore, the next-generation transcriptome data from 42 echinoderm specimens from 24 orders, 37 families have been collected to establish a web-based application (http://echinotol.org) to identify orthologs suitable for phylogenetic analysis in assembling the echinoderm tree of life .
In sea cucumber, transcriptomic sequencing has been performed together with genomic and proteomic analyses in Apostichopus japonicus, to facilitate the molecular underpinnings of visceral regeneration . RNA-seq has been applied independently to analyze the DEGs between normal and regenerating radial nerve cord in Holothuria glaberrima, and revealed the key roles of extracellular matrix (ECM) remodeling and ECM-cell interactions in regeneration . In addition, the applications of transcriptomic analyses in sea cucumber include identification of long noncoding RNA species , illustration of the mechanisms of aestivation , abnormal development , and body wall pigmentation . Moreover, immune-related genes in the A. japonicus coelomocytes under Vibrio splendidus challenge have been identified, which are clustered into the immune pathways of endocytosis, lysosome, chemokine, and MAPK and ERBB signaling . However, there is still limited research reports on the species and response of the cytokines secreted by coelomocytes of sea cucumber under challenge of immune stimuli.
In this study, high-throughput transcriptomic sequencing and bioinformatics analysis were performed on the primary coelomocytes isolated from the sea cucumber H. leucospilota and treated with three different immune stimuli including lipopolysaccharide (LPS), polyinosinic-polycytidylic acid [poly (I:C)] and heat-inactived Vibrio harveyi for 24 h. Based on the DEGs between different groups, the immune-related pathways were screened out and the responded cytokines were identified. This study provides evidences for the potential roles of cytokines in the innate immunity of sea cucumber.
Illumina draft reads and sequence assembly
Primary coelomocytes isolated from the sea cucumber H. leucospilota were respectively challenged with LPS, poly (I:C) and heat-inactivated V. harveyi for 24 h (Fig. 1a), and twelve cDNA libraries were constructed to perform Illumina sequencing. After assembly and redundancy removal, a transcriptome with twelve RNA-seq libraries in total of 6.69 GB with 73,472 identified Unigenes was obtained and submitted to GenBank under the BioProject accession No. PRJNA559679. The total length, average length, N50, and GC content of the Unigenes were 47,163,631 bp, 641 bp, 1015 bp and 39.54%, respectively. The number of transcripts and Unigenes decreased with increasing of length, and the majority of them were concentrated in 200–3000 bp (Fig. 1b). The Unigene sequences were annotated to seven functional databases, and 20,926 (28.48%) of them were significantly matched to at least one of the databases (Table 1), in which 19,156 (26.07%) to NR, 3990 (5.43%) to NT, 14803 (20.15%) to SwissProt, 13,615 (18.53%) to KOG, 15277 (20.79%) to KEGG, 7097 (9.66%) to GO and 15,624 (21.27%) to InterPro. The annotation results for five databases were further showed in a Venn diagram: 1584, 11, 88, 21 and 528 genes were independently annotated into the NR, KOG, KEGG, SwissProt and InterPro databases, respectively, and the intersection set of these five databases was 11,103 (Fig. 1c). For the NR database, the annotated Unigenes were majorly matched to Strongylocentrotus purpuratus (5552, 28.98%), Acanthaster planci (6117, 31.93%), Saccoglossus kowalevskii (883, 4.61%), Branchiostoma belcheri (457, 2.39%) and other species (6147, 32.09%, Fig. 1d).
Analysis of DEGs and validation of RNA-Seq data by RT-qPCR
The transcriptomic data obtained from the control group (CT), LPS-treatment group (LPS), Poly (I:C)-treatment group (PIC) and V. harveyi-treatment group (VH) were analyzed comparatively. The results showed that the co-expressed DEGs for the three groups were 1180, and the uniquely-expressed genes in the LPS-vs-CT, PIC-vs-CT and VH-vs-CT groups were 3846, 3869 and 2279, respectively (Fig. 2a). The significantly DEGs were acquired by comparing of the gene expression between the LPS, PIC and VH groups and the CT group with the following criteria: P ≤ 0.01, |log 2-fold-change| ≥ 1 and false discovery rate (FDR) ≤ 0.05. Finally, 7074 DEGs in the LPS-vs-CT comparison (666 upregulated and 6408 downregulated), 7737 DEGs in the PIC-vs-CT comparison (355 upregulated and 7382 downregulated), and 5481 DEGs in the VH-vs-CT comparison (387 upregulated and 5094 downregulated) were obtained (Fig. 2b).
To validate the gene expression results of the RNA-seq data, 13 significant DEGs from the LPS-vs-CT comparison with the criteria of P < 0.01, FDR < 0.05 and fold change of at least 4 were selected for RT-qPCR validation. The volcano plot for the CT-vs-LSP comparison was shown in Fig. 2c. The expression levels of the selected 13 genes in the LPS challenge group were normalized to the control group. As anticipated, the RT-qPCR data showed a positive linear relationship with the RNA-Seq data (Fig. 2d), and there was no statistically significant difference between the two datasets (P > 0.05). This result suggested that the RNA-Seq was a positively related reference for expression profiling study on the whole, and the assembly quality of the sequences was desirable.
Functional classification and enrichment analysis of GO terms and KEGG pathways for co-expressed DEGs
The co-expressed DEGs after the three different immune challenges may represent the essential genes for immune defense in sea cucumber. GO (gene ontology) is a major bioinformatics initiative to unify the classification of genes and gene product attributes across all species. Among the total 1180 co-expressed DEGs, 796 were annotated with GO terms, in which 321 were annotated as “biological process”, 279 were annotated as “cellular components” and 196 were annotated as “molecular functions” (Fig. 3a).
KEGG (Kyoto Encyclopedia of Genes and Genomes) is a reference database dealing with genomes, biological pathways, diseases, drugs, and chemical substances. In this case, 1003 co-expressed DEGs were annotated in KEGG, and the number of DEGs in cellular processes, environmental information processing, genetic information processing, human diseases, metabolism and organismal systems categories were 123, 113, 20, 344, 161 and 242, respectively (Fig. 3b).
Gene set enrichment analysis (GSEA, also called functional enrichment analysis) is a method to identify classes of genes or proteins that are overrepresented in a large set of genes or proteins and . The co-expressed DEGs in this study were annotated into 273 KEGG pathways. The top 10 and top 20 significantly enriched pathways were shown in Table 2 and Fig. 3c, respectively. The pathways of Staphylococcus aureus infection, complement and coagulation cascades, and tubercuosis were the three most significantly enriched pathways, with a total of 15, 16 and 24 co-expressed DEGs annotation, respectively (Table 2). In addition, the top 10 immune-related KEGG pathways with the highest number of co-expressed genes were shown in Table 3. Within them, the most significant immune-related KEGG pathways were complement and coagulation cascades, which contained 84 annotated genes with 16 co-expressed DEGs. These results laid a foundation for discovering the immune-related genes and developing the immune responding mechanisms of sea cucumber. The functional classification and enrichment analysis of GO terms and KEGG pathways for the DEGs were further performed for the challenges of different immune stimuli (Fig. 3d). In this case, the totals of 2591, 3260 and 2585 Unigenes respectively for the comparisons of the LPS-vs-CT, PIC-vs-CT and VH-vs-CT were annotated in GO terms. Among three different immune challenge groups, the Poly (I:C) treatment group had the most DEGs annotated into the GO terms, indicating that more genes were involved in the administration of Poly (I:C) in the sea cucumber coelomocytes when compared to other two immune stimuli.
The most significant enrichment pathway of DEGs in different immune challenges
A hypergeometric test was used for enrichment analysis of all signal pathways in the KEGG database. Compared to the control group, the number of DEGs annotated to specific pathways for the LPS, Poly (I:C) and V. harveyi treatment groups were 3415, 4261 and 3120, respectively. The top 10 most significantly enriched pathways for each challenge were shown in Table 4, in which the common pathways for the three treatments groups were the ABC transporter pathway, ubiquitin-mediated proteolysis pathway and inositol phosphate metabolism pathway. The glucagon signaling pathway was a common pathway for the Poly (I:C) and V. harveyi treatment groups, and other pathways were specific for each treatment. Analysis of the DEGs enriched KEGG pathways provides an effective basis for studying the immune defense process, biological function and metabolic pathways in the sea cucumber coelomocytes.
Analysis of immune-related pathway in different immune challenges
Compared to the control group, the DEGs of the LPS, Poly (I: C) and V. harveyi treatment groups were categorized into 308, 316 and 305 annotated KEGG pathways, respectively. Accordingly, 18, 22 and 23 significantly enriched immune-related pathways for LPS, Poly (I: C) and V. harveyi treatments, respectively, were identified according to the immune system of KEGG pathways. The most significantly enriched immune-related pathways for LPS, PIC and V. harveyi treatments were the Th1 and Th2 cell differentiation signaling pathways, the intestinal immune network for IgA production pathway, and the Fc gamma R-mediated phagocytosis pathway, respectively. Other top 10 significantly enriched immune-related pathways for the three immune challenges are shown in Table 5. Base on the analysis of different regulatory pathways of the DEGs in KEGG, the mechanism of cellular immune response of the sea cucumber coelomocytes can be understood more directly.
Identification of cytokines and their expression analysis after different challenges
Cytokines are a broad and loose category of small proteins (~ 5–20 kDa) that are important for cell signaling, especially the immune signaling. Twenty-one cytokines were selected according to the NR database annotation of the transcriptomic data. The identified cytokines belong to four cytokine families, namely, the B-cell lymphokine (BCL/CLL), erythroid differentiation-related factor 1-like (EPRF1), interleukin 17-like (IL-17) and thrombospondin-like (TSP/TPO) families. In our transcriptome for the sea cucumber H. leucospilota, the BCL/CLL family included CCL11A, BCL3-X3, BCL10, CLL7A, CLL9-X3a, CLL9-X3b, CLL9-X3c, and CLL9-X3d; the EPRF1 family included EDRF1a, EDRF1b, EDRF1c, and EDRF1d; the IL-17 family included IL-17, IL-17-2, IL-17B, and IL-17C/E; and the TSP/TPO family included TPO-I-7A1–3, TPO-I-7B1–3, TSP-1a-c, and TSP-4. Among these cytokine families, the expression levels of the four genes in the interleukin-17 family after LPS challenge were significantly higher than those of other cytokines (Fig. 4), indicating that IL-17 was an important family of the cytokines in the immune response of sea cucumber, and its effective mechanism needs to be further investigated.
Phylogenetic and structural domain analysis of the selected cytokines
A phylogenetic tree was constructed using a maximum-likelihood (ML) method under MEGA7.0 with the deduced amino acid sequences from the selected cytokines, including BCL/CLL, EPRF1, IL-17 and the TSP/TPO family (Fig. 5a). The result showed that the cytokines were clustered into the corresponding branches.
The structural domains for some of the selected cytokines were predicted and shown in Fig. 5b. Specially, the members of TSP/TPO family contain repeated type 1 thrombospondin domains. The BCL/CLL family cytokines contain a conserved N-terminal domain that is found in the BCL7 family. The predicted IL-17, IL-17-2, IL-17B, and IL-17C/E in the IL-17 family commonly contain a cysteine knot fold domain.
Transcript expression of cytokines in the coelomocytes with dosage-increased LPS treatment
The selected cytokines were applied to study their expressing upon LPS treatments of increased concentrations of 10.0 μg/ml, 20.0 μg/ml and 50.0 μg/ml. As the result, the mRNA levels of IL-17 cytokines family were significantly upregulated after 10 μg/ml LPS challenge, and their expression levels generally showed a trend of first increasing, then decreasing, and finally stabilizing with the increasing of LPS concentration (Fig. 6). In the BCL/CLL family, with the increase of LPS concentration, BCL10, CLL7A, CLL9-X3b and CLL9-X3c showed a dose-dependent increasing expression pattern, similar to the expression patterns of TPO-I-7A2 and TPO-I-7A2 in the TSP/TPO family. In contrast, CCL11A, BCL3-X3, CLL9-X3a and CLL9-X3d were expressed in a dose-dependent descending pattern, similar to the expression trends of TPO-I-7B1–3 and TSP-4 in the TSP/TPO family. In the EPRF1 family, EDRF1a, EDRF1b, EDRF1c, and EDRF1d showed similar parallel expression patterns, which first increased, then decreased, and finally stabilized with increasing of LPS concentrations (Fig. 6).
Based on Illumina HiSeq2500 sequencing, a transcriptome for the sea cucumber primary coelomocytes under challenges of LPS, Poly (I: C) and V. harveyi was obtained in this study. Without referenced genomic data, the assembled Unigene sequences were annotated using the databases of nucleotides and genomes of other species. There are 28.98 and 31.91% of the unigenes that annotated to the NR database were matched to the sequences from S. purpuratus and A. planci, respectively. One of the possible reason is that the first two reference genomes databases were supplied in these two species [27, 28]. However, other two studies for the genome of the sea cucumber A. japonicas have been published recently [19, 29], which can provide more conveniences for annotating the genes in the transcriptome of other sea cucumber speices in the future.
A Unigene in the transcriptomic data is declared differentially expressed if a difference observed in read counts between the experimental and control conditions is statistically significant. Thirteen Unigenes were selected from the DEGs for RT-qPCR verification (Fig. 2c & d). The results showed that the trends of gene expression by RT-qPCR were consistent with those by RNA-seq analysis, verifying the credibility of the RNA- seq results.
The co-expressed DEGs in all three groups with different immune challenge [LPS, Poly (I: C) and V. harveyi] may speculated to be the essential genes of the immune defense in sea cucumber, while the DEGs unique to a kind of immune stimuli may specifically respond to the corresponding stimuli. LPS is the main component of the wall of gram-negative bacteria, as a typical endotoxin, can bind to the CD14/TLR4/MD2 complexes in many cell types, especially monocytes, dendritic cells, macrophages and B cells, which can promote the secretion of pro-inflammatory cytokines, nitric oxide and eicosenoic acid, leading to a series of immune responses in the host . In the LPS-treatment group, the Th1 and Th2 cell differentiation pathway was the most significantly enriched immune-related pathway. Poly (I:C), structurally similar to double-stranded RNA , is an immune activator that is used to mimic viral infection, and it is known to interact with TLR3, which is expressed on the cell membrane of B cells, macrophages and dendritic cells. In the Poly (I:C) treatment group, the intestinal immune network for IgA production pathway was the most significantly enriched immune-related pathway. V. harveyi is a gram-negative marine bacterium that is pathogenic to commercial aquaculture species, including shrimp, fish, and sea cucumber . In the V. harveyi treatment group, the Fc gamma R-mediated phagocytosis pathway was the most significantly enriched immune-related pathway. Annotation of the DEGs into the immune-related pathways provides a basis for the research of immune molecular mechanism of echinoderms.
For analysis of DEGs in the coelomocytes with different immune challenges, the results showed that compared with the control group, more down-regulated DEGs than up-regulated DEGs were presented in the groups challenged with three different immune stimuli (Fig. 2b), indicating that the expression for a majority of genes were suppressed by immune challenge. The top 10 significantly enriched KEGG pathways (Tables 2 and 4) and immune-related KEGG pathways (Tables 3 and 5) for common and specific immune challenged DEGs were further screened out. Based on the transcriptomic data in this study, the cDNAs of H. leucospilota Fas-associated death domain (FADD) have been cloned with functional characterization  It has been reported that there are many humoral immune factors in the coelomic fluid of echinoderms that can recognize and attack the invaded pathogens, including different types of lectin, lysozyme, hemolysin, hydrolase, phenol oxidase, superoxide dismutase, nitric oxide synthase and complement-like factor [6, 9, 34,35,36]. Similarly, a number of humoral immune factors in the corresponding immune-related pathways (Table 5) were identified in the current study.
The top 10 KEGG pathways for the most significant DEG enrichment after different immune challenges also included amino acid and carbohydrate metabolic pathways, such as alanine, aspartate and glutamate metabolism pathway, amino sugar and nucleotide sugar metabolism pathway, and glycosaminoglycan biosynthesis-heparan sulfate/heparin pathway (Table 4). A similar study performed in in the sea cucumber A. japonicas has reported that the arginine metabolic pathway is related to the pathogenic challenge with a dose-dependent manner . Combined with our current study, it is speculated that the metabolism of some amino acids is related to the host immune response in echinoderm, and the regulatory network between immunity and metabolism is more complicated than known now.
Based on the DEGs identified in the transcriptome of sea cucumber coelomocytes challenged with immune stimuli, twenty-one candidate genes of immune cytokines were selected for analyzing the cytokine response of sea cucumber against invaded pathogens. The selected cytokines belong to four cytokine families, namely, the TSP/TPO, BCL/CLL, EPRF1 and IL-17 families (Fig. 5a). Among these gene families, the expression levels of 4 genes in the IL-17 family after LPS challenge were significantly higher than those of other cytokines. By further confirmation with RT- qPCR of the expression levels of 4 IL-17 s generally showed a trend of first increasing, then decreasing, and finally stabilizing with incubation of LPS with increasing concentrations (Fig. 6). IL-17 is a kind of potent proinflammatory cytokine produced by activated memory T cells in mammals . The IL-17 family (of which there are 6 known members, termed IL-17A to IL-17F) is thought to represent a distinct signaling system that appears to highly conserved across vertebrate evolution . IL-17 family members play an active role in inflammatory diseases, autoimmune diseases and cancer . Based on the identification of differentially expressed cytokine genes in the sea cucumber coelomocytes under immune challenges, this study may provide an evidence that many kinds of cytokines interact with each other, and play a complex and essential role in the innate immunity of sea cucumber.
At present, more and more studies on the immune system of the sea cucumber have been developed at both the molecular and cellular levels, which were based on the basal information for immune-related genes that provided by the transcriptomic studies . However, the specific immune defense mechanism in the effector coelomocytes is still unclear. In this study, high-throughput transcriptome sequencing and bioinformatics analysis were performed on coelomocytes treated with three different immune stimuli (LPS, Poly (I:C) and Vibrio harveyi). Analysis of the results could elucidate the molecular mechanism for the immune response in the sea cucumber coelomocytes to the pathogenic/immune challenge, laying a foundation for the future identification of immune-related genes and characterizations of immune-responded mechanisms of echinoderm.
In this study, comparative transcriptome analysis of DEGs among LPS (10 μg/ml), Poly (I:C) (10 μg/ml) and heat-inactivated Vibrio harveyi (107 cells/ml) immune stimuli of sea cucumber H. leucospilota coelomocytes revealed the DEGs and immune-related pathways that are crucial for research of molecular mechanisms related to the immune response in echinoderms. Twenty-one cytokine candidate DEGs were identified, which belong to 4 cytokine families, namely, BCL/CLL, EPRF1, IL-17 and TSP/TPO. Among them IL-17 family cytokines were significantly upregulated after 10 μg/ml LPS challenge for 24 h in response to LPS dose-increased treatment (0, 10, 20 and 50 μg/ml), which provide new insights for the echinoderm cytokine response during immune challenge.
Animals used in this research were obtained from commercial sea cucumber catches, therefore approval from any ethics committee or institutional review board was not necessary. A total of nine healthy H. leucospilota weighing 190–210 g were obtained from the Daya Bay in Guangdong province, China, and maintained in a seawater aquarium with aerated and filtrated seawater (30‰ salinity) at 32 °C for 1 week before the experiments were performed.
Isolation, primary culture and static incubation of coelomocytes
Given that coelomocytes are considered as the immune effector cells in echinoderms , the coelomocytes from H. leucospilota were applied to explore the response of sea cucumber after challenge of immune stimuli. Isolation of primary cells and static incubation of coelomocytes were performed as previously reported with some modifications [40, 41]. Briefly, three sea cucumbers were washed with sterile DEPC water three times before dissection on ice. The coelomic fluids were sieved through a 150-mesh cell sieve to remove large tissue debris, mixed with the prechilled anticoagulant solution (20 mM EGTA, 480 mM NaCl, 19 mM KCl and 68 mM Tris, pH 7.6) in a 1:1 volume ratio. The cell suspension was filtered through a 100-μm nylon mesh, and the cells were collected by centrifugation at 500×g at 4 °C for 10 min. The harvested cells were washed twice with 30 ml isotonic buffer (1 mM EGTA, 530 mM NaCl and 10 mM Tris, pH 7.6) and resuspended in 10 ml Leiboviz’s L-15 cell culture medium (Gibco BRL, USA; containing 390 mM NaCl, 100 U/mL penicillin and 100 μg/ml streptomycin, pH 7.6) and filtered through a 40-μm nylon mesh to remove the cell clusters. The cell yield was counted with a 0.10-mm hematocytometer, and only the cell preparations with greater than 95% viability, assessed by trypan blue exclusion assay, were used in subsequent experiments.
The cell suspensions were diluted to 0.1 × 106 cells/ml and seeded onto two 24-well culture plates (1 ml/well) that were precoated with ploy D-lysine. The isolated coelomocytes were incubated at 28 °C for another 18 h for recovery. On the following day, 1 ml of immune stimulus prepared in Leiboviz’s L-15 medium was gently overlaid onto the cells after removal of the old medium. The immune stimuli used in this study included: LPS (10 μg/ml), Poly (I:C) (10 μg/ml) and heat-inactivated Vibrio harveyi (107 cells/ml), and fresh Leiboviz’s L-15 medium was used as a control. For RNA sequencing, the cells were harvested at 24 h after incubation of tested substrates. In this case, cells from four wells were pooled together as a biological duplicate, and three duplicates for each group were sequenced and analyzed.
RNA extraction, library construction and sequencing
Total RNA from each sample was extracted using TRIzol Reagent (Invitrogen, USA). The quality and concentration of total RNA were determined with NanoDrop (Thermo Scientific, USA), and the RNA integrity value (RIN) was checked using the RNA 6000 Pico LabChip on an Agilent 2100 Bianalyzer (Agilent, USA). The RNA-seq library was constructed in six steps: 1) ployA-tailed mRNAs were enriched by oligo (dT) selection; 2) the obtained mRNAs were fragmented and reverse transcribed into double-stranded cDNA (dscDNA) by N6 random primer; 3) the cDNA ends were repaired by 3′-adenylation and adaptor ligation; 4) the ligation products were amplified by PCR; 5) the PCR products were heat-denatured and the single-stranded DNAs were cyclized; and 6) sequencings were performed on an Illumina HiSeq2500 (Illumina, USA).
De novo assembly and Unigene functional annotation
The raw sequencing data were first filtered by using SOAPnuke (v1.5.2) software under the parameters of -l 15 -q 0.2 -n 0.1 to remove the reads in which unknown bases (N) are more than 5%. The low quality reads, which were defined by the base qualities for more than 20% of them were less than 15, were then removed. After reads filtering, de novo assembly was performed using Trinity2.0.6 with clean reads . Tgicl2.0.6 was used on cluster transcripts to remove redundancy and get Unigenes . After assembly, the Unigenes functional annotation were performed with seven functional databases, namely NR, NT, GO, KOG, KEGG, SwissProt and InterPro. Specifically, We use Blastn 2.2.23 , We use Blastx 2.2.23  or Diamond 0.8.31  align Unigenes to NT, NR, KOG, KEGG and SwissProt database to do the annotation, use Blast2GO 2.5.0  with NR annotation to do the GO annotation, and use InterProScan5 5.11–51.0  to do the InterPro annotation, under default parameter.
Quantification of gene expression level and identification of DEGs
Clean reads of each sample were mapped to the Unigenes with Bowtie2 2.2.5 , and the gene expression levels were calculated with RSEM 1.2.12 (RNA-Seq by Expectation-Maximization) to obtain the FPKM (Fragments Per Kilobase of exon model per Million mapped reads) value . The DEGseq algorithm was used to detect the differences of gene expression between different groups. To improve the accuracy of DEGs, genes with over two-fold changes (log2 ratio ≥ 1) and q-value < 0.001 were considered as significantly DEGs.
Screening and analysis of immune-related pathways
Based on Wallenius non-central hyper-geometric distribution , GO enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseqR packages 1.10.0, which can be adjusted for gene length bias in DEGs. KOBAS (KEGG Orthology-Based AnnotationSystem) software was used to test the statistical enrichment of DEGs in KEGG pathways. The FDR correction was performed for p-value, and pathways with FDR < 0.01 were regarded as significantly enriched pathways. The top 10 immune-related pathways significantly enriched in the groups treated with different stimuli were selected by the immune-related annotation of the KEGG pathways.
Validation of differentially expressed genes by RT-qPCR
Significant DEGs (P < 0.05, FDR < 0.05, fold change ≥8) were selected for RT-qPCR analysis to validate the transcriptomic data, and corresponding primers were designed by based on the Unigenes sequences (Table 1). The primary coelomocytes isolated from another three H. leucospilotas were then challenge by LPS (10 μg/ml) for 24 h. Total RNA from each sample was extracted using TRIzol Reagent and reverse transcribed into the first cDNA using Honorll 1st strand cDNA Synthesis SuperMix for qPCR kit (Novogene, China). Finally, quantitative real-time PCRs were performed using the Unique Aptamer qPCR SYBR probe kit (Novogene). In this case, the qPCRs were performed in triplicates for cells cultured in three individual wells.
Identification and bioinformatics analysis of cytokines responding to immune challenge
According to the NR annotation of the transcriptomic data from the coelomocytes challenged by LPS, Poly (I:C) and V. harveyi, Twenty-one candidate cytokine genes were selected. The log2-fold changes of their expression levels under different challenges were presented as heatmap by the Pheatmap package using R3.5 software [51, 52].
The CDSs of the candidate cytokine genes were translated to amino acid sequences to construct a maximum-likelihood phylogenetic tree with 1000 bootstrap replicates by using MEGA 7. The structural domains of the cytokines were predicted by using the SMART program.
Measurement of cytokines mRNA levels in coelomocytes after LPS challenge
Primers for the selected cytokine genes were designed by based on their Unigenes sequences (Table 6). The primary coelomocytes isolated from the other three H. leucospilotas were obtained after treatment of LPS in gradient increased concentrations (0, 10, 20 and 50 μg/ml) for 24 h. Total RNA extraction, first-strand cDNA synthesis and RT-qPCR were performed as described above. In this case, β-actin was used as a reference gene, and the RT-qPCRs were performed in triplicates for cells cultured in three individual wells.
Availability of data and materials
RNA-seq data from this article have been deposited in NCBI under the accession number of STUDY: PRJNA559679, SAMPLE: D2CT (SAMN12551561), EXPERIMENT: C1_1 (SRX6693868), RUN: D2CT1_1.fq.gz (SRR9945372), SAMPLE: D2LPS (SAMN12551562), EXPERIMENT: L2_2 (SRX6693869), RUN: D2LPS1_1.fq.gz (SRR9945371), SAMPLE: D2PIC (SAMN12551563), EXPERIMENT: P3_3 (SRX6693870), RUN: D2PIC1_1.fq.gz (SRR9945370), SAMPLE: D2VH (SAMN12551564), EXPERIMENT: V4_4 (SRX6693871), RUN: D2VH1_1.fq.gz (SRR9945369).
Analysis of variance
Differentially expressed genes
Erythroid differentiation-related factor 1
Fragments per kilobase of exon model per million mapped reads
Gene set enrichment analysis
Kyoto encyclopedia of genes and genomes
Clusters of orthologous groups for eukaryotic complete genomes
Mitogen-activated protein kinase
- Poly (I:C):
RNA integrity value
RNA-seq by expectation-maximization
Real time quantitative polymerase chain reaction
Bonham K, Held EE. Ecological observations on the sea cucumbers holothuria atra and H. leucospilota at Rongelap Atoll, Marshall Islands. University of Hawaii Press; 1963. p. 305-14.
Huang W, Huo D, Yu ZH, Ren CH, Jiang X, Luo P, Chen T, Hu CQ. Spawning, larval development and juvenile growth of the tropical sea cucumber Holothuria leucospilota. Aquaculture. 2018;488:22–9.
Blair JE, Hedges SB. Molecular phylogeny and divergence times of deuterostome animals. Mol Biol Evol. 2005;22(11):2275–84.
Qian J, Ren CH, Chen T, Xia JJ, Yu ZH, Gao Y, Hu CQ. Identification and functional characterization of first alpha-2-macroglobulin in sea cucumbers. Aquac Res. 2017;48(5):2278–90.
Qian J, Ren CH, Xia JJ, Chen T, Yu ZH, Hu CQ. Discovery, structural characterization and functional analysis of alpha-2-macroglobulin, a novel immune-related molecule from Holothuria atra. Gene. 2016;585(2):205–15.
Chen T, Ren CH, Li WH, Jiang X, Xia JJ, Wong NK, Hu CQ. Calmodulin of the tropical sea cucumber: gene structure, inducible expression and contribution to nitric oxide production and pathogen clearance during immune response. Fish Shellfish Immunol. 2015;45(2):231–8.
Yan AF, Ren CH, Chen T, Jiang X, Sun HY, Huo D, Hu CQ, Wen J. The first tropical sea cucumber caspase-8 from Holothuria leucospilota: molecular characterization, involvement of apoptosis and inducible expression by immune challenge. Fish Shellfish Immunol. 2018;72:124–31.
Yan AF, Ren CH, Chen T, Huo D, Jiang X, Sun HY, Hu CQ. A novel caspase-6 from sea cucumber Holothuria leucospilota: molecular characterization, expression analysis and apoptosis detection. Fish Shellfish Immunol. 2018;80:232–40.
Xue Z, Li H, Wang XL, Li X, Liu Y, Sun J, Liu CJ. A review of the immune molecules in the sea cucumber. Fish Shellfish Immunol. 2015;44(1):1–11.
Ramirez-Gomez F, Aponte-Rivera F, Mendez-Castaner L, Garcia-Arraras JE. Changes in holothurian coelomocyte populations following immune stimulation with different molecular patterns. Fish Shellfish Immunol. 2010;29(2):175–85.
Schillaci D, Cusimano MG, Cunsolo V, Saletti R, Russo D, Vazzana M, Vitale M, Arizza V. Immune mediators of sea-cucumber Holothuria tubulosa (Echinodermata) as source of novel antimicrobial and anti-staphylococcal biofilm agents. AMB Express. 2013;3:35.
Wang TT, Sun YX, Jin LJ, Xu YP, Wang L, Ren TJ, Wang KL. Enhancement of non-specific immune response in sea cucumber (Apostichopus japonicus) by Astragalus membranaceus and its polysaccharides. Fish Shellfish Immunol. 2009;27(6):757–62.
Mehdi Jazayeri S, Melgarejo Munoz LM, Mauricio Romero H. RNA-SEQ: a glance at technologies and methodologies. Acta Biolo Colomb. 2015;20(2):23–35.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.
Wygoda JA, Yang Y, Byrne M, Wray GA. Transcriptomic analysis of the highly derived radial body plan of a sea urchin. Genome Biol Evol. 2014;6(4):964–73.
Israel JW, Martik ML, Byrne M, Raff EC, Raff RA, McClay DR, Wray GA. Comparative developmental transcriptomics reveals rewiring of a highly conserved gene regulatory network during a major life history switch in the sea urchin Genus Heliocidaris. PLoS Biol. 2016;14(3):e1002391.
Fuess LE, Eisenlord ME, Closek CJ, Tracy AM, Mauntz R, Gignoux-Wolfsohn S, Moritsch MM, Yoshioka R, Burge CA, Harvell CD, et al. Up in arms: immune and nervous system response to sea star wasting disease. PLoS One. 2015;10(7):e0133053.
Janies DA, Witter Z, Linchangco GV, Foltz DW, Miller AK, Kerr AM, Jay J, Reid RW, Wray GA. EchinoDB, an application for comparative transcriptomics of deeply-sampled clades of echinoderms. BMC Bioinformatics. 2016;17(1):48.
Zhang XJ, Sun LN, Yuan JB, Sun YM, Gao Y, Zhang LB, Li SH, Dai H, Hamel JF, Liu CZ, et al. The sea cucumber genome provides insights into morphological evolution and visceral regeneration. PLoS Biol. 2017;15(10):31.
Mashanov VS, Zueva OR, Garcia-Arraras JE. Transcriptomic changes during regeneration of the central nervous system in an echinoderm. BMC Genomics. 2014;15:357.
Mu C, Wang R, Li T, Li Y, Tian M, Jiao W, Huang X, Zhang L, Hu X, Wang S, et al. Long non-coding RNAs (lncRNAs) of sea cucumber: large-scale prediction, expression profiling, non-coding network construction, and lncRNA-microRNA-gene interaction analysis of lncRNAs in Apostichopus japonicus and Holothuria glaberrima during LPS challenge and radial organ complex regeneration. Mar Biotechnol. 2016;18(4):485–99.
Zhao Y, Yang H, Storey KB, Chen M. RNA-seq dependent transcriptional analysis unveils gene expression profile in the intestine of sea cucumber Apostichopus japonicus during aestivation. Comp Biochem Physiol Part D Genomics Proteomics. 2014;10:30–43.
Li Y, Kikuchi M, Li X, Gao Q, Xiong Z, Ren Y, Zhao R, Mao B, Kondo M, Irie N, et al. Weighted gene co-expression network analysis reveals potential genes involved in early metamorphosis process in sea cucumber Apostichopus japonicus. Biochem Biophys Res Commun. 2018;495(1):1395–402.
Xing L, Sun L, Liu S, Li X, Zhang L, Yang H. De novo assembly and comparative transcriptome analyses of purple and green morphs of Apostichopus japonicus during body wall pigmentation process. Comp Biochem Physiol Part D Genomics Proteomics. 2018;28:151–61.
Gao QLM, Wang Y, et al. Transcriptome analysis and discovery of genes involved in immune pathways from coelomocytes of sea cucumber (Apostichopus japonicus) after Vibrio splendidus challenge. Int J Mol Sci. 2015;16(7):16347–77.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
Sea Urchin Genome Sequencing C, Sodergren E, Weinstock GM, Davidson EH, Cameron RA, Gibbs RA, Angerer RC, Angerer LM, Arnone MI, Burgess DR, et al. The genome of the sea urchin Strongylocentrotus purpuratus. Science. 2006;314(5801):941–52.
Hall MR, Kocot KM, Baughman KW, Fernandez-Valverde SL, Gauthier MEA, Hatleberg WL, Krishnan A, McDougall C, Motti CA, Shoguchi E, et al. The crown-of-thorns starfish genome as a guide for biocontrol of this coral reef pest. Nature. 2017;544(7649):231–+.
Li YL, Wang RJ, Xun XG, Wang J, Bao LS, Thimmappa R, Ding J, Jiang JW, Zhang LH, Li TQ, et al. Sea cucumber genome provides insights into saponin biosynthesis and aestivation regulation. Cell Discov. 2018;4:29.
Glinski Z, Jarosz J. Immune phenomena in echinoderms. Arch Immunol Ther Exp. 2000;48(3):189–93.
Fortier ME, Kent S, Ashdown H, Poole S, Boksa P, Luheshi GN. The viral mimic, polyinosinic : polycytidylic acid, induces fever in rats via an interleukin-1-dependent mechanism. Am J Phys Regul Integr Comp Phys. 2004;287(4):R759–66.
Li Y, Xu X-L, Zhao D, Pan L-N, Huang C-W, Guo L-J, Lu Q, Wang J. TLR3 ligand poly IC attenuates reactive Astrogliosis and improves recovery of rats after focal cerebral ischemia. CNS Neurosci Ther. 2015;21(11):905–13.
Zhao L, Ren C, Chen T, Sun H, Wu X, Jiang X, Huang W. The first cloned sea cucumber FADD from Holothuria leucospilota: molecular characterization, inducible expression and involvement of apoptosis. Fish Shellfish Immunol. 2019;89:548–54.
Shao Y, Che Z, Xing R, Wang Z, Zhang W, Zhao X, Jin C, Li C. Divergent immune roles of two fucolectin isoforms in Apostichopus japonicus. Dev Comp Immunol. 2018;89:1–6.
Wang H, Xue Z, Liu ZQ, Wang WL, Wang FF, Wang Y, Wang LL, Song LS. A novel C-type lectin from the sea cucumber Apostichopus japonicus (AjCTL-2) with preferential binding of D-galactose. Fish Shellfish Immunol. 2018;79:218–27.
Jiang J, Zhou Z, Dong Y, Zhao Z, Sun H, Wang B, Jiang B, Chen Z, Gao S. Comparative expression analysis of immune-related factors in the sea cucumber Apostichopus japonicus. Fish Shellfish Immunol. 2018;72:342–7.
Austin B, Zhang XH. Vibrio harveyi: a significant pathogen of marine vertebrates and invertebrates. Lett Appl Microbiol. 2006;43(2):119–24.
Aggarwal S, Gurney AL. IL-17: prototype member of an emerging cytokine family. J Leukoc Biol. 2002;71(1):1–8.
Kolls JK, Linden A. Interleukin-17 family members and inflammation. Immunity. 2004;21(4):467–76.
Ren C, Chen T, Jiang X, Wang Y, Hu C. Identification and functional characterization of a novel ferritin subunit from the tropical sea cucumber, Stichopus monotuberculatus. Fish Shellfish Immunol. 2014;38(1):265–74.
Ren C, Chen T, Jiang X, Luo X, Wang Y, Hu C. The first echinoderm gamma-interferon-inducible lysosomal thiol reductase (GILT) identified from sea cucumber (Stichopus monotuberculatus). Fish Shellfish Immunol. 2015;42(1):41–9.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–U130.
Pertea G, Huang XQ, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B, et al. TIGR gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics. 2003;19(5):651–2.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12(1):59–60.
Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.
Quevillon E, Silventoinen V, Pillai S, Harte N, Mulder N, Apweiler R, Lopez R. InterProScan: protein domains identifier. Nucleic Acids Res. 2005;33:W116–20.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–U354.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Fog A. Calculation methods for Wallenius’ noncentral hypergeometric distribution. Commun Stat Simul Comput. 2008;37(2):258–73.
Marwick BBC, Mullen L. Packaging data analytical work reproducibly using R (and friends). PeerJ Preprints. 2018;6:e3192v2.
Ihaka R. R: past and future history, vol. 30; 1998.
We are thankful to Mr. Zezhong, Zhang (Sino-Danish College, University of Chinese Academy of Sciences) and Dr. Hongmei, Li (Institution of South China Sea Ecology and Environmental Engineering, Chinese Academy of Sciences) for help with submitting data.
This study was supported by the National Key R & D Program of China (2018YFD0901605), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA13020205), the Science & Technology Promoting Projects for Oceanic & Fishery in Guangdong Province (SDYY-2018-01) and the Guangdong Province Program (2017B030314052, 2018A030313857, 2015A030310120, A2015230). The National Natural Science Foundation of China (41676162), the Science and Technology Service Network Initiative of Chinese Academy of Sciences (KFJ-STS-ZDTP-055). The Science and Technology Research Programs of Guangdong Province (2018A030313857). The funding bodies had no role in the design of the study, the collection, analysis, and interpretation of data, and in writing the manuscript and the decision to submit the work for publication.
Ethics approval and consent to participate
Animals used in this research were obtained from commercial sea cucumber catches, therefore approval from any ethics committee or institutional review board was not necessary. No ethics approval was required for the collection or experimentation of the animal used in this study, the Holothuria leucospilota.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Wu, X., Chen, T., Huo, D. et al. Transcriptomic analysis of sea cucumber (Holothuria leucospilota) coelomocytes revealed the echinoderm cytokine response during immune challenge. BMC Genomics 21, 306 (2020). https://doi.org/10.1186/s12864-020-6698-6
- Sea cucumber
- Differentailly expressed genes
- Immune response