Transcriptome Analysis Reveals the Tolerance Mechanism of Mantis Shrimp (Oratosquilla oratoria) under a Lipopolysaccharide Challenge.

Lipopolysaccharide (LPS), a major cell wall component of Gram-negative bacteria, is considered to lead to some disease development in commercial crustaceans. However, mantis shrimps Oratosquilla oratoria (Crustacea: Stomatopoda) have a strong vitality and ability to resist disease. To study the tolerance mechanism of mantis shrimp, transcriptome analyses were conducted in hepatopancreas of O. oratoria under LPS challenge investigation. Totally, 84 547 044 clean reads were obtained from transcriptomes (43 159 230 in OP (control), 41 387 814 in OL (treatment), respectively). Unigenes, the longest transcript of each gene, with a total length of 68 318 880 bp and the total number of 100 978 were obtained. 8369 (8.28%) of unigenes were successfully annotated in all databases and 54 888 (54.35%) were annotated in at least one database. Finally, 1012 differentially expressed genes (DEGs) including 439 and 573 showed significantly upregulated and downregulated were determined between OL and OP, respectively. Moreover, those DEGs only expressed in OL or OP accounted for 8.99%. The functional classification based on GO and KEGG indicated that the common enrichment categories for the DEGs are "amino sugar metabolic" and "cellular homeostasis" and that the progress of nutrient metabolic and homeostasis in cells is important in facing variable environmental conditions. Protein-protein interaction analysis elucidated proteins, β-actin (ACTB_G1), T-complex protein subunits (TCPs), heat shock proteins (HSPs), hydroxysteroid dehydrogenase-like protein 2 (HSDL2), kinesin family member 5 (KIF5), methylglutaconyl-CoA hydratase (AUH), and myosin heavy chain (MYH) may play key roles in response to an LPS challenge. This study laid a foundation to further investigate the possible adaptation way that O. oratoria survives in a bacterial challenge.

Introduction

Diseases in aquaculture are increasing in both frequency and intensity, and for crustaceans, a recent report mentions that up to 40% of total shrimp production was lost annually mainly due to disease impacts.[1] Lipopolysaccharide (LPS), a major cell wall component of Gram-negative bacteria, was considered to lead to some disease development in crustaceans.[2,3] Some crustacean diseases, such as early mortality syndrome (EMS), milky hemolymph disease (MHD), and necrotizing hepatopancreatitis (NHP) were caused by bacterial agents.[1] These diseases also represented the toxic effects of bacterial pathogens. It is important to understand the toxicity mechanism of bacterial pathogens in crustaceans, especially those with high economic values.

LPS, a key component of pathogenic Vibrio, has been widely utilized in studies on the toxic effects, the toxicity mechanism in crustaceans,[1,4−8] and the response to toxicity stress.[9−12] For example, pathogenic Vibrio has been shown to colonize the gastrointestinal tract of shrimps and produce a toxin (a highly positively charged protein manufactured by the EMS strains of Vibrio parahaemolyticus) that can cause tissue destruction and dysfunction of the hepatopancreas.[1,13] Meanwhile, it has been indicated that haemocyte apoptosis would eliminate damaged or weak cells and contribute to haemocyte renewal in response to toxicity from pathogens and that it may be a defense strategy against pathogens in tiger shrimp Penaeus monodon.[11] Recently, the development of high-throughput RNA-sequencing studies had made some efforts based on “Big Data” from transcriptome to reveal the response mechanism in response to various stimuli including LPS.[14−17] It is known that mainly important immune-related pathways were involved in response to LPS including MAPK, NF-κB, and Jak-STAT signaling pathways.[18] However, some other pathways may also refer to an LPS challenge.[19,20] Moreover, several important proteins in MAPK, NF-κB, and Jak-STAT signaling pathways had also been elucidated to play a role in response to LPS in the Eriocheir sinensis transcriptome.[3] Transcriptome has been applied in various crustaceans, such as Litopenaeus vannamei, Sinopotamon henanense, Triops longicaudatus, Scylla paramamosain, and Exopalamon carincauda.[14,21−23] Thus, transcriptome analysis may be a powerful tool to reveal some mechanisms in response to various stimuli from the environment.

Mantis shrimp (Oratosquilla oratoria) is one of the commercially valuable aquatic species in the Indo-West Pacific region, and they have a strong vitality and ability to resist bacterial infection and heavy metals.[24,25] However, less information about the tolerance mechanism to pathogenic bacteria was available until now on O. oratoria. Therefore, to reveal the antibacterial mechanism on this species, this study investigated mantis shrimp’s responses to an LPS challenge at the transcriptome level.

Materials and Methods

Sample Collection and Treatment

Twenty individuals O. oratoria (length: 11.8 ± 0.5 cm) were collected from the aquatic market in October 2016 in Yancheng, Jiangsu Province, China. Then, they were acclimatized at 20° centigrade for 72 h prior to the natural environment. Twenty individuals were divided into two groups (ten control groups (OP) and ten treatment groups (OL)). All treatment groups were injected with 5 mg/kg of LPS (L-2654, Sigma) in 15 μL of PBS, and the control groups were injected with 15 μL PBS. The hepatopancreas in two groups was collected after 24 h, frozen in liquid nitrogen, and stored at −80 °C for later experiments.

RNA Isolation and Library Construction

Total RNA was extracted using TRIzol reagent (Invitrogen Corp., Carlsbad, CA), and its concentration was measured using the Qubit RNA Assay Kit in a Qubit 2.0 Flurometer (Life Technologies, CA) according to the manufacturer’s instructions. A total amount of 1.5 μg RNA per hepatopancreas sample was used in this study. The integrity of RNA was assessed using the RNA Nano 6000 Assay Kit of an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA). Libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB) following manufacturer’s instructions, and index codes were added to attribute sequences to each sample.

Sequencing and De Novo Transcriptome Assembly

The libraries were sequenced on an Illumina HiSeq. 2500 platform. Raw data (raw reads) in fastq format were first processed using in-house perl scripts. Next, clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N, and low-quality reads. Transcriptome assembly was performed using the Trinity program.[26] During assembly, sequences were clustered by Corset to remove redundant sequences, and these sequences were defined as transcripts.[27] Finally, unigenes were produced from the transcript sequences according to the longest transcript for each gene.

Functional Unigene Annotation and Classification

Unigenes were annotated based on the following databases: Nr/Nt (NCBI nonredundant protein/nucleotide sequences) (http://www.ncbi.nlm.nih.gov/), Pfam (Protein family) (https://pfam.sanger.ac.uk/), Swiss-Prot (A manually annotated and reviewed protein sequence database) (http://www.ebi.ac.uk/uniprot), KOG/COG (Clusters of Orthologous Groups of proteins) (https://www.ncbi.nlm.nih.gov/COG/), KO (KEGG Ortholog database) (http://www.genome.jp/kegg/), and GO (Gene Ontology) (http://www.geneontology.org/). To obtain significant annotations, alignments in these databases use 10–5 as the cut-off E-value. The further classification of functional unigenes were conducted mainly by Gene Ontology (GO) terms using Blast2GO and KEGG.[28,29]

Determination of Differentially Expressed Genes (DEGs) and Their Enrichment

The determination of DEGs between the two groups was performed using the DESeq R package.[30] The resulting p-values were adjusted using the Benjamini and Hochberg method to control the false discovery rate. Genes with an adjusted p-value <0.05 as identified by DESeq were assigned as differentially expressed. A Gene Ontology (GO) enrichment analysis of the DEGs was performed using the GOseq R package based on Wallenius’ noncentral hypergeometric distribution, which can adjust for gene length bias in DEGs. We used KOBAS software to test for statistical enrichment of the DEGs in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.[31]

Analysis of Protein–Protein Interactions (PPIs)

The sequences of the DEGs were blasted against the genome of a related species (the protein interactions of which can be found in the STRING database [http://string-db.org/]) to predict PPIs. The PPI network of the DEGs was visualized in Cytoscape.[32]

Results

Sequencing and De Novo Transcriptome Assembly

A total of 84 547 044 clean reads were obtained in the two transcriptomes (43 159 230 in OP and 41 387 814 in OL); 59.33% (25 604 604) and 65.56% (27 132 468) of the clean reads in OP and OL, respectively, were successfully mapped onto the reference (i.e., the sequences that assembled clean reads). Further de novo transcriptome assembly based on hierarchical clustering using Corset showed that the total length and number of transcripts were 108 261 683 bp and 244 699, respectively. Moreover, 100 978 unigenes were obtained with a total length of 68 318 880 bp. The length of an assembled transcriptome (e.g., the maximum, average, and minimum) indicates the reliability of the assembly, especially the N50 value, which serves as a relative standard to determine the quality of the transcriptome assembly; these characteristics are provided in Table .[33] In this study, the N50 value was determined to be 489 and 801 bp for transcripts and unigenes, respectively; these values are not small compared with similar studies.[34,35] The most abundant transcripts were clustered in a group of <301 bp in length, the most abundant unigenes were clustered in a group of 501–1000 bp in length, and 2963 unigenes were clustered in a group of >2000 bp in length (Figure ). The assembly results also indicate that the length distribution pattern and mean length of the unigenes are consistent with those in previous transcriptome studies using the Illumina HiSeq platform.[22,36]

Table 1

Summary of the De Novo Assembly of the Transcriptome Profiles of O. oratoria

sample	total nucleotides	no.	max length (bp)	ave length (bp)	N50 (bp)
transcripts	108 261 683	244 699	25 305	301	489
unigene	68 318 880	100 978	25 305	515	801

Figure 1

Length distribution of the unigene and transcript of O. oratoria.

Unigene Annotation and Functional Classification

Annotation of the 100 978 unigenes in seven databases (National Center for Biotechnology Information [NCBI] nonredundant protein [Nr], NCBI nucleotide sequences [Nt], KEGG Orthology [KO], GO, Swiss-Prot, protein family [Pfam], and EuKaryotic Orthologous Groups [KOG]) showed that 8369 (8.28%) unigenes were successfully annotated in all databases and 54 888 (54.35%) were annotated in at least one database. Thus, more than 40% of the unigenes could not be annotated. Most of the unigenes were annotated in GO (37 686, 37.32%) followed by Nr (36 963, 36.61%) and Swiss-Prot (36 948, 36.59%) (Table ).

Table 2

Summary of Annotation of 100 978 Unigenes in O. oratoria

database	number of unigenes	percentage (%)
annotated in Nr	36 963	36.61
annotated in Nt	35 325	34.98
annotated in KO	18 583	18.4
annotated in SwissProt	36 948	36.59
annotated in Pfam	36 753	36.39
annotated in GO	37 686	37.32
annotated in KOG	19 488	19.29
annotated in all databases	8369	8.28
annotated in at least one database	54 888	54.35
total unigenes	100 978	100

Functional classification of the unigenes was achieved by GO and KEGG analyses. In the GO analysis, the annotated unigenes (37 686) were classified into three categories: “cell component”, “biological process”, and “molecular function”. These categories included 165 662, 509 444, and 79 556 unigenes, respectively (Figure ). In the KEGG analysis, most of the unigenes (18 583) analyzed were assigned to “metabolism” (7444) compared with “cellular processes” (3354), “environmental information processing” (2645), “genetic information processing” (4897), and “organismal systems” (6170) (Figure ). Genes like NfKb,[37] a member of the RAS oncogene family (RAP1A), and activating transcription factor 1 (ATF),[38] which have been reported to respond to lipopolysaccharide (LPS), were also annotated in this study. Annotation and functional classification of the unigenes laid the foundation for further functional enrichment analysis of DEGs.

Figure 2

Gene Ontology (GO) annotations of the assembled unigenes: 3 categories and 56 sub-categories.

Figure 3

Cluster of KEGG Ortholog database annotations of assembled unigenes.

DEGs and Functional Enrichment

A total of 1012 DEGs were selected according to the q-value (adjusted p-value) and |log2 (fold change)| value (>1). Of the 1012 DEGs, 439 and 573 were significantly upregulated and downregulated, respectively, between OL and OP (Figure ). Moreover, the DEGs expressed only in OL or OP accounted for a total of 8.99% (Table S1); these DEGs play important roles in adaption to an LPS challenge. In this study, the functional enrichment of DEGs by GO analysis first showed that the DEGs were primarily significantly enriched in root categories with the terms “biological process” and “molecular function” with an adjusted p-value of <0.05 (Figure ). Descendent GO terms included the “amino sugar metabolic process”, “chitin metabolic process”, “cellular metal ion homeostasis”, and “chitin binding”. Functional enrichment based on KEGG pathway analysis revealed that these DEGs could be assigned to 201 pathways, and the number of DEGs ranged from 1 to 17. Of these pathways, 25 had a p-value <0.05 (Table S2), indicating that according to Hui, 25 pathways might be associated with the response to an LPS challenge.[15] The top 20 pathways were selected, and 6 pathways (“protein digestion and absorption”, “ECM-receptor interaction”, “neuroactive ligand-receptor interaction”, “focal adhesion”, “amebiasis”, and “amino sugar and nucleotide sugar metabolism”) had adjusted p-values of <0.05 (Figure ). Thus, genes involved in these pathways may be important in response to an LPS challenge. Finally, six upregulated and eight downregulated pathways were significantly enriched (adjusted p-value <0.05).

Figure 4

Volcano plot of differentially expressed genes between OP and OL.

Figure 5

GO terms of the DEGs significantly enriched in O. oratoria hepatopancrea.

Figure 6

Statistics of KEGG pathway enrichment analysis between OP and OL.

PPIs of DEGs

PPI analysis allowed us to further characterize the DEGs. DEGs with a higher “degree” are more likely to play a key role in response to an LPS challenge, as a higher degree suggests a more important function.[39] In this study, higher degree DEGs included those encoding β-actin (ACTB_G1), T-complex protein subunits (CCT6, CCT8, and CCT5, also known as TCPs), heat shock proteins (HSPA5 and HSPA1_8), hydroxysteroid dehydrogenase-like protein 2 (HSDL2), and kinesin family member 5 (KIF5) (Table and Figure S1). ACTB_G1 was assigned to “focal adhesion” with an adjusted p-value of <0.05.

Table 3

Profiles of DEGs with Higher Degrees in the KEGG Pathway

unigenes	description	KO name
Cluster-42114.8042	β-actin	ACTB_G1
Cluster-42114.9886	T-complex protein 1 subunit zeta	CCT6
Cluster-42114.7141	T-complex protein subunit theta	CCT8
Cluster-42114.9040	heat shock protein cognate 3	HSPA5
Cluster-42114.11094	heat shock protein	HSPA1_8
Cluster-42114.8274	heat shock protein 70	HSPA1_8
Cluster-10940.0	T-complex protein subunit γ	CCT3
Cluster-42114.2893	hydroxysteroid dehydrogenase-like protein 2	HSDL2
Cluster-30781.1	kinesin family member 5	KIF5
Cluster-17651.1	methylglutaconyl-CoA hydratase, mitochondrial precursor	AUH
Cluster-42114.2274	T-complex protein 1 subunit epsilon	CCT5
Cluster-42114.9439	myosin heavy chain	MYH

Discussion

Cell matrix adhesion plays essential roles in important biological processes, including cell motility, proliferation, differentiation, gene expression regulation, and survival (ko04510). Although ACTB_G1 is involved in basic housekeeping functions, varying expression levels of ACTB_G1 in response to different stimuli, including exposure to madrigal (downregulation), herpes simplex virus type 1 (downregulation), hormones and serum (upregulation), have been reported.[40] However, whether upregulation of ACTB_G1 occurs in response to an LPS challenge remains unclear. T-complex proteins (CCT6, CCT8, and CCT5) and heat shock proteins, which act as chaperonins, function not only in proper protein folding but are also involved in immune responses.[41] Although its function is generally unknown, HSDL2 may function as a regulatory factor in fatty acid metabolism.[42,43] In human glioblastoma cell lines, HSDL2 knockdown significantly inhibited cell proliferation, resulting in cancer cell apoptosis.[43] In O. oratoria, upregulation of HSDL2 might be a negative response to an LPS challenge. Similar to HSDL2, upregulation of methylglutaconyl-CoA hydratase (AUH) may also be a negative response to an LPS challenge. AUH is present in mitochondria, with an enriched distribution in organs with higher energy demands; however, both the knockdown and overexpression of AUH have been linked to a decreased abundance of respiratory complexes.[44] The kinesin KIF5, upregulated in “endocytosis” with a p-value <0.05, also caught our attention. However, the functional significance of members of the kinesin family is only partially understood.[45] To date, the primary function of kinesins is in the axonal transport of membrane proteins and KIF5 has been shown to regulate endosomal membrane traffic to mediate a variety of physiological functions.[45−47] Therefore, KIF5 may be indirectly involved in the response to an LPS challenge.[48] In this study, myosin heavy chain (MHC) with five direct edges in PPI analysis was also investigated. Unlike vertebrates, in which the function of MHC has been well studied, only limited information is available regarding invertebrate MHC. For example, one type of human MHC has been shown to suppress the heat-induced death of airway smooth muscle cells.[49] In crustaceans, studies have focused on the characterization of muscle, since myosin is the major protein found in skeletal muscle.[48] Fortunately, research into disease pathology based on natural history data and experimental interventions has indicated that changes in MHC expression might confer some resistance to damage to compensate for dystrophic dysfunction.[50] Upregulated MHC in O. oratoria in response to an LPS challenge might be associated with this protective mechanism.

The response of the mantis shrimp to an LPS challenge at the transcriptome level indicates the following: First, the common enrichment categories of the DEGs were “amino sugar metabolic” and “cellular homeostasis”, which complied with the objective rule that metabolic processes and homeostasis in cells are important in facing variable environmental conditions. Second, the pathways “protein digestion and absorption”, “ECM-receptor interaction”, “neuroactive ligand-receptor interaction”, “focal adhesion”, and “amebiasis” with adjusted p-values of <0.05 should be investigated in future studies. Third, our PPI analysis revealed that ACTB_G1, TCPs, HSPs, HSDL2, KIF5, AUH, and MYH might play key roles in response to an LPS challenge. Real-time RT-PCR analysis showed that the DEGs annotated as ACTB_G1, CCT6, CCT8, HSPA5, HSDL2, and KIF5 were upregulated in the treatment group compared to the control (Figure ).

Figure 7

Results of representative DEGs between RT-PCR and RNA-seq.

Nuclear factor (NF)-κB is a member of a family of transcription factors that function as dimers and regulate genes involved in immunity, inflammation, and cell survival. The canonical NF-κB pathway is induced by tumor necrosis factor-α, interleukin-1, or byproducts of bacterial and viral infections.[19,51−53] The “NF-kappa B signaling pathway” in this study was enriched for two upregulated unigenes (Cluster-42114.8893 and Cluster-42114.14556), which were annotated as COX2 and NFKB1.[20] When we combined the GO and KEGG enrichment results, the common enrichment categories of the DEGs were amino sugar metabolic and cellular homeostasis, which were partly consistent with a proteomic study of the sea cucumber Apostichopus japonicas (“amino sugar metabolic”).[20] In immune cells, metabolic processes and cellular homeostasis are important; these cells patrol the body via the bloodstream and migrate into different tissues where they face variable environmental conditions, different substances, and changes in oxygen levels.[54] In the present study, the common enrichment categories comply with the objective rule that O. oratoria responded to an LPS challenge.