Single-Cell Sequencing Reveals Types of Hepatopancreatic Cells and Haemocytes in Black Tiger Shrimp (Penaeus monodon) and Their Molecular Responses to Ammonia Stress.

The cell types and developmental trajectories of shrimp cells based on the transcriptional level have not been established, and gene expression profile and function at the single-cell level is unclear. We aimed to use scRNA-seq to construct a single-cell resolution transcriptional map of hepatopancreas and haemocytes in shrimp to analyse the molecular mechanisms of the immune response to ammonia nitrogen stress. In the present study, seven cell clusters were successfully identified in each of the two tissues (haemocytes, Hem1-7; hepatopancreas, Hep1-7) based on specifically-expressed marker genes. The developmental starting points of haemocytes and hepatopancreatic cells were Hem2 and Hep1, respectively. We propose that Hem2 has oligopotent potential as the initiation site for haemocyte development and that Hem4 and Hem5, located at the end of development, are the most mature immune cell types in haemocytes. Hep5 and Hep6 were the developing terminal cells of hepatopancreas. The antioxidant system and proPO system of shrimp were activated under ammonia nitrogen stress. A large number of DEGs were involved in oxidative stress, detoxification metabolism, and immune defence. In particular, important response genes such as AMPs, proPO, and GST were not only marker genes for identifying cell groups but also played an important role in shrimp cell differentiation and functional plasticity. By successfully applying 10× Genomics based scRNA-seq to the study of shrimp, the single-cell transcriptional profiles of hepatopancreatic cells and haemocytes of shrimp innate immune responses under ammonia stress were constructed for the first time. This atlas of invertebrate hepatopancreatic cells and haemocytes at single-cell resolution identifies molecular events that underpin shrimp innate immune system responses to stress.

Introduction

Crustaceans are a morphologically and physiologically diverse animal group and are one of the most successful animal taxa worldwide. In decapod crustaceans, the hepatopancreas is the largest organ of the digestive gland and is the main metabolic organ, which is responsible for the absorption and metabolism of nutrients, digestion and storage, synthesis of digestive enzymes and proteins, and detoxification of xenobiotics (1). Although the hepatopancreas has been investigated using histological, physiological, and molecular biological techniques there is no universally accepted structural and functional conceptualization (2).

Shrimp are the main focus of crustacean aquaculture. Haemocytes play an important role in the immune defence system of shrimp; they execute cellular immune processes and synthesise and release several immune factors to provide the material basis for humoral immunity. Shrimp haemocytes are divided into three morphological types (hyaline, semi-granular, and granular cells) based on their morphology and physiology (3, 4). The establishment of shrimp cell lines has been a challenge and the function of individual cell types is unknown. The full functional diversity of shrimp haemocytes and their developmental trajectories has not been established, and the extent of functional equivalence in morphologically-similar haemocytes is unclear. However, characteristics of these cells are largely determined by cells’ molecular signature or gene expression patterns. Therefore, in order to further understand shrimp cell diversity and function, it is necessary to identify the cell types and their marker genes.

Ammonia nitrogen (ammonia-N) is a serious pollutant that has a negative impact on aquaculture systems. Moreover, ammonia-N accumulation can induce a series of changes in histology, metabolic characteristics, and immune system of organisms, leading to neurological dysfunction, impaired respiration, decreased fecundity, and animal diseases and death, which can cause large economic losses to farmers (5). Efforts have been made to decipher the molecular mechanisms of the detrimental effects of ammonia on aquatic animals, such as shrimp and crabs (6, 7). Haemocytes and hepatopancreas are the key tissues that respond to ammonia stress in shrimp, and their functions are related to metabolism, detoxification, and immune defence (5, 8). However, the responding mechanism of haemocytes and the hepatopancreas in shrimp has been studied mainly at a whole tissue or bulk RNA sequencing levels. Due to the limited understanding of cellular diversity, the transient state of cells and their dynamic gene expression patterns are often buried in the overall RNA sequencing data. Hence, it is important to construct a single-cell molecular atlas of important tissues of shrimp, characterizing the molecular signatures of all dynamic states of cell types under steady state and stress conditions.

Single-cell mRNA sequencing (scRNA-seq) is a powerful tool for studying the state and function of cells. Recently, scRNA-seq has been used to reveal cellular characteristics and functions in invertebrates, dramatically changing the understanding of the cellular diversity of multiple invertebrates, such as mosquitoes, Drosophila, silkworms, and scallops (9–12). In particular, as one of the mainstream platforms for scRNA-seq technology, 10x Genomics Chromium system has the advantages of high throughput and high cell capture rate, and is widely used in the research of cell heterogeneity, cell classification, construction of cell maps, and discovery of new cell types. Here, we performed scRNA sequencing of haemocytes and hepatopancreas of P. monodon, one of the most common commercial crustacean species, under steady state and ammonia stress, to comprehensively distinguish mature cell types from their transient intermediate states. We focused on developing an atlas of diverse cell populations in shrimp (haemocytes and hepatopancreas) to provide a valuable resource for the gene expression profiles of the various cell types and their states in shrimp.

Material and Methods

Animals and Ammonia-N Exposure Experiments

Healthy P. monodon (8.2 ± 1.0 g) were obtained from the South China Sea Fisheries Research Institute in Shenzhen (Guangdong, China). 120 shrimp were randomly selected and equally divided into two groups (control, CG; ammonia-stressed, AG) with two replicates each. A replicate was one 500-L fibreglass tank containing 30 shrimp. In the AG group, NH4Cl solution was added until tanks reached the predetermined ammonia concentration. The ammonia concentration was set to ∼58.15 mg/L (LD50 for 96 h of ammonia exposure). Water salinity was 30 ppt; temperature, 27 ± 1°C; pH, 7.8–8.1; and dissolved oxygen, ≥ 6.5 mg/L. After 12 h, three shrimp per group were separately dissected to obtain hepatopancreas tissues and haemolymph.

Preparation of Hepatopancreatic Cells and Haemocytes

Fresh hepatopancreas tissue were immediately used to prepare single-cell suspensions. Briefly, the collected hepatopancreas tissue was washed with DPBS (PBS without calcium and magnesium) and carefully chopped into small pieces. Enzymatic digestion was performed with collagenase II (2 mg/mL) and dispase (0.2 mg/mL) at 37°C for 15min. The digested solution was transferred to a centrifuge tube after passing through a 70 μm filter and a 40 μm stainless nylon mesh (Greiner Bio-OneGmbH, Germany). The cell suspension was centrifuged at 400 g for 4 min, and the supernatant was discarded slowly. The cell layer was resuspended in DPBS containg 0.04% BSA after washing with PBS at 400 g for 4 min twice.

To avoid agglutination, anticoagulant EDTA solution was used in haemocytes separation. The needle and syringe were wet with EDTA prior to collection of haemolymph. The collected haemolymph was mixed with EDTA solution in a 1:1 volume. Then the collected haemolymph was immediately transferred in a 50 ml centrifuge tube filled with 20 ml DPBS containing 0.04% BSA, and centrifuged at 300 g for 5 min to collect the haemocytes. The collected haemocytes were washed twice with PBS and centrifuged at 400 g for 4 min and resuspended in DPBS containing 0.04% BSA. Cell viability was measured using a Trypan Blue staining kit (Sangon Biotech). The cells of each sample were counted using a cell counting plate.

Single-Cell RNA Sequencing

The library was synthesised as per the instructions for Chromium Next GEM Single Cell 3′ Reagent Kits v3.1 (10× Genomics). In brief, the cells from three shrimp of each group were mixed to form one sample and adjusted to 1000 cells/μL. Cell suspensions were loaded on the Chromium Controller instrument to generate single-cell Gel Beads-In-Emulsions (GEMs). In GEMs, individual cells were barcode-labelled and GEM-reverse transcriptions were performed. Then, the cDNA libraries were amplified using primers of R1 and P5 arms after reverse transcription. Sequence-ready libraries were constructed according to the Chromium Single Cell 3′ Library v3 Kit (10× Genomics) instructions. Finally, four barcoded sequencing libraries (two for hepatopancreatic cells and two for haemocytes) were sequenced on the Illumina 10× Genomics Chromium platform.

Quality Control (QC) and Data Normalisation

The raw FASTQ data were processed using CellRanger software (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/5.0/what-is-cell-ranger), which performs splicing-aware alignment of reads to P. monodon genome (13) by calling STAR (2.7.2a) (https://github.com/alexdobin/STAR) with default parameters. To retain high-quality single-cell RNA-seq data, genes expressed in fewer than three cells, low-quality cells with gene numbers less than 200, cells with gene numbers higher than 2500, and cells with more than 0.05% mitochondrial genes, were removed in the initial QC screen by the CellRanger. Then, doublets in raw data were identified and removed by DoubletFinder (https://github.com/chris-mcginnis-ucsf/DoubletFinder) within a range of 7.5%. After removing unwanted cells from the dataset, the “Normalization” function of the Seurat software was used to normalise the expression of data (https://satijalab.org/seurat/pbmc3k_tutorial.html).

Cell Clustering and Visualisation

To cluster the cells, principal component analysis (PCA) was performed using the normalised expression amount. The most significant ∼12–15 principal components (PCs) were selected from the PCA results for subsequent cluster analysis, and the single-cell sub-group classification results were further visualised by Uniform Manifold Approximation and Projection (UMAP) (https://satijalab.org/seurat/articles/get_started.html) in a two-dimensional space.

Differentially Expressed Genes (DEGs) Analysis

Differential gene expression analysis on different cell populations was performed using Seurat’s bimod likelihood ratio statistical test (https://satijalab.org/seurat/articles/de_vignette.html). The DEGs were determined by the following criteria: p-value ≤ 0.05, average log2FC ≥ 0.25, and the expressed gene detection in more than 25% of cells.

Visualisation of Marker Genes

To illustrate the gene expression patterns, heatmaps and dot plots for all the cell clusters were constructed using marker genes. The marker genes were the upregulated genes in a single cluster compared to their expression in all other cells, which were screened using a non-parametric Wilcoxon rank sum test with a threshold of adjusted p-value < 0.05. The top 10 upregulated genes in each cluster were used for plotting heatmaps.

Functional Enrichment Analyses

Functional enrichment analyses were performed for the DEGs and marker genes. For GO enrichment analysis, DEGs were first mapped to the InterProScan database (https://www.ebi.ac.uk/interpro/search/sequence/) to obtain a list of genes with a certain GO function in each term, and significantly enriched GO terms were identified by a hypergeometric test with a threshold of p-value ≤ 0.05. For KEGG enrichment analysis, a Clusterprofiler package (v4.0.5) (https://github.com/YuLab-SMU/clusterProfiler) was applied to find pathways significantly enriched in DEGs compared to the background of the entire transcript set with a threshold of p-value ≤ 0.05.

Pseudotime Inference

To explore the temporal resolution of cell differentiation states, single cell trajectory was performed using Monocle 2 (https://github.com/cole-trapnell-lab/monocle-release). Then, a pseudotime ordering of cells is inferred from the DEGs. We assumed that the raw UMI counts were distributed according to a negative binomial distribution with fixed variance expression family to model the raw UMI count data (14, 15).

Results

Global Transcriptional Profiling and Cell Type Identification of Adult Shrimp Cells

We obtained 12,712 haemocytes and 21,226 hepatopancreatic cells for subsequent scRNA-seq analysis after rigorous QC ( ). We detected a median value of 5072 UMIs and 1222 genes per cell in haemocytes, and 1733 UMIs and 645 genes per cell in hepatopancreatic cells. The two types of integrated single-cell transcriptomes were dimensionally reduced through UMAP. Based on highly variable genes, the haemocyte’s were found to contain seven main cell clusters (Hem1-Hem7; ) with a proportion range of 1.35% (Hem7) to 38.02% (Hem1). Seven hepatopancreas cell clusters (Hep1-Hep7) were identified and their proportions ranged from 0.83% (Hep7) to 43.23% (Hep1) ( ). The details of marker genes of each cluster of haemocytes and hepatopancreatic cells are listed in .

Figure 1

ScRNA-seq of P. monodon haemocytes and hepatopancreas reveals their cellular subpopulations. (A) Overall workflow for cell sorting and single-cell data analyses. (B) UMAP plot of cell clusters from shrimp haemocytes. (C) UMAP plot of cell clusters from shrimp hepatopancreas. (D) Heatmap of the top 10 marker genes in haemocytes. (E) Heatmap of the top 10 marker genes in hepatopancreas. (F) Violin plots of the significantly marker genes in haemocytes. (G) Violin plots of the significantly marker genes in hepatopancreas.

Heatmaps were constructed to display the top 10 most variably-expressed genes for each cluster in haemocytes ( ) and hepatopancreas ( ) respectively, and the most specific genes from each cell cluster were visualized using violin plots in both haemocytes ( ) and hepatopancreas ( ). The results showed that each of the seven clusters had specific biomarkers which highlighted dynamic gene expression changes occurring throughout adult shrimp haemocytes and hepatopancreas.

Differential Cell Number Between AG and CG Groups and the Marker Gene Enrichment Atlas

UMAP analysis displayed cell cluster composition and their changes in haemocytes and hepatopancreatic cells ( ). The comparison of cell number ratios between CG and AG conditions showed that the clusters of Hem1 and Hem6 were significantly increased, while clusters of Hem4 and Hem7 were significantly decreased in the AG compared to that in the CG condition ( ). For hepatopancreatic cells, the numbers of Hep1, Hep3, Hep5, and Hep6 cells were significantly increased, while the number of Hep7 cells was significantly decreased while the number of Hep7 cells was significantly decreased after ammonia-N stress ( ).

Figure 2

The changes of cell clusters under ammonia-stress and enrichment of marker genes in P. monodon haemocytes and hepatopancreas (A) UMAP of haemocytes in the control (CG) and ammonia-treated groups (AG). The bubble chart represents the changes in the cell ratio of each cell cluster in haemocytes. (B) UMAP of hepatopancreas in the AG and CG. The bubble chart represents the changes in the cell ratio of each cell cluster in hepatopancreas. (C) GO and KEGG enrichment analysis of marker genes in haemocytes. (D) GO and KEGG enrichment analysis of marker genes in hepatopancreas.

Further, GO and KEGG enrichment analyses of the marker genes with significant expression differences in each cluster ( ) revealed that the marker genes of different clusters in haemocytes were mainly related to immune phagocytosis, homeostasis, and detoxification. For example, clusters of Hem1 and Hem2 were mainly related to peptide synthesis, amino acid metabolism, and carbon metabolism. The marker genes of the Hem3/4/5 clusters were mainly enriched in homeostasis, peptide metabolism, amino acid activation, Ras GTPase binding, actin cytoskeleton organisation, lipid metabolism, and regulation of response to stimulus. The involved KEGG pathways were glutathione metabolism, endocytosis, phagosome, and autophagy ( ). The marker genes of the Hem6/7 clusters were mainly enriched in cell cycle and death-related events, such as DNA replication, DNA repair, and apoptosis ( ). Hem7 is a cell cluster related to phagocytosis and apoptosis, and maker genes were mainly enriched in molecular pathways such as endocytosis, apoptotic process and cytoskeleton organization.

In terms of cell clusters derived from hepatopancreatic tissues ( ), the clusters of Hep1/2/3/6 cells were mainly enriched in metalloendopeptidase activity, drug metabolism, ATP metabolism, aromatic compound catabolism, and apoptosis ( ). In general, the biological functions of the Hep1/2/3/6 cell clusters were mainly associated with immune phagocytosis and detoxification; and Hep6 was associated with apoptosis. Parts of the marker genes of Hep5 were related to organic acid metabolism and inorganic anion transport. In addition, the marker genes of the Hep4/7 clusters were enriched in GO terms of DNA packaging, ATP metabolism, and regulation of cell death. The Hep4 and Hep7 clusters were associated with chromosome replication events, metabolism, and energy regulation ( ).

Shrimp Developmental Trajectory in the Haemolymph and Hepatopancreas

Furthermore, RNA velocity was used to predict the dynamic changes in cell development in shrimp haemocytes and hepatopancreas ( ). The haemocytes appeared to develop along two differentiation pathways, among which the developmental trajectory of Path I was Hem 2-Hem 3-Hem 4 (end point) and the trajectory of Path II was Hem2-Hem1-Hem5 (end point; ). We applied monocle2 to perform pseudotime analysis of Hem1-Hem5. According to the change in trajectory, the cells underwent three stages: the starting point of the branch (pre-branch) and the other two branches (cell fate1 and cell fate2; ). The genes enriched in pathways such as ribosome synthesis, were highly expressed at the pre-branch stage (P1; ). The expression levels of most enriched genes related to metabolism and oxidation pathways were significantly increased at cell fate2 (P2; ). However, ECM receptor interaction, autophagy and other pathways were significantly increased at cell fate1 but significantly decreased at cell fate2 (P3; ). In addition, the expression patterns in haemocytes of several important genes were notable ( ). For example, glutathione peroxidase, GSH-Px (LOC119594415) was highly expressed in cell fate1. AMPs (LOC119594995 and LOC119594996) were highly expressed in cell fate2. The expression levels of anti-lipopolysaccharide factor, ALF (LOC119576449) and prophenoloxidase, proPO (LOC119586224) were higher in cell fate2 ( ).

Figure 3

The pseudotemporal trajectories identify transcriptional dynamics of shrimp haemocytes and hepatopancreatic cells. (A) UMAP diagram of the RNA velocity of haemocytes. (B) UMAP diagram of the RNA velocity of hepatopancreatic cells. (C) Pseudotime analysis of cells in Path I and Path II development directions. (D) Pseudotime analysis of cells in Hep1 and Hep3 developmental directions. (E) The pseudotime branch heatmap of haemocytes and KEGG enrichment. (F) The pseudotime branch heatmap of hepatopancreatic cells and KEGG enrichment. (G) Monocle differentiation tree of Hep1/5/6 in hepatopancreas. (H) The branch heatmap of Hep1/5/6 cells and KEGG enrichment. (I) Jitter diagram of branch development of key genes in haemocytes. (J) Jitter diagram of key genes in the developmental direction of hepatopancreatic cells.

For hepatopancreatic cells, RNA velocity analysis showed that the differentiation directions were mainly Path I (Hep1-Hep3), II (Hep1-Hep5-Hep6) and III (Hep2). Two developmental pathways of Path I and Path II were selected for pseudotime analysis ( ). During the development process of Path I, the highly expressed genes in the start stage were mainly enriched in pathways such as lysosome, starch and sucrose metabolism, and neuroactive ligand–receptor interaction (P2; ), and were mainly enriched in TGF signal transduction, endocytosis and RNA transport at the end stage (P1; ). The pseudotime analysis results of Path II showed three main stages based on the trajectory plot ( ). Lysosome and starch and sucrose metabolism pathways were highly expressed at the pre-branch stage ( ). Ribosome and endocytosis pathways were significantly decreased at cell fate1 but significantly increased at cell fate2 (P2; ). In addition, we also noted the expression patterns in the hepatopancreas of several important genes ( ). For example, the perlucin-like protein (LOC119574611) was the most highly expressed in the pre-branch state and decreased with development. The expression level of vesicle-associated membrane protein 3 (LOC119583416) was increased at the cell fate1.

The Differential Molecular Profiles of Haemocytes Between AG and CG Groups

To gain a deeper understanding of the potential mechanisms of the blood lymphocyte response to ammonia stress, the DEGs of the three focused clusters of Hem3, Hem4, and Hem5 between the AG and CG groups were analysed ( ). The DEGs of Hem3 were enriched in drug metabolism and the mTOR signalling pathway. The DEGs of Hem4 and Hem5 were enriched in oxidative phosphorylation; the DEGs of Hem3 and Hem5 were enriched in mitophagy; and the DEGs of these three cell clusters were all enriched in the ribosome, phagosome, and endocytosis pathways ( ). The intergroup differences of upregulated and downregulated genes in the three clusters are shown in a Venn diagram ( ). Of the 243 co-upregulated genes, 60.91% were functionally annotated ( Down). Of the 114 co-downregulated genes, 78.95% were functionally annotated ( Down). Violin plots were used to demonstrate the characteristic genes upregulated and downregulated in response to ammonia stress in the Hem3, Hem4, and Hem5 clusters ( ).

Figure 4

Enrichment analysis of differentially expressed genes in P. monodon haemocytes under ammonia-stress. (A) Differential gene volcanic map display of Hem3/4/5 clusters in haemocytes. (B) KEGG enrichment analysis of differential genes in Hem3/4/5 clusters of haemocytes. (C) The Venn diagram of coupregulated genes and its gene function in haemocytes. (D) The Venn diagram of codownregulated genes and its gene function in haemocytes. (E) Violin diagram showing functional genes in haemocytes. (*p < 0.05 and **P < 0.01).

The Differential Molecular Profiles of AG and CG Groups of Hepatopancreas

For the differential analysis of shrimp hepatopancreatic cells, we selected the clusters that developed from Hep1 to Hep3, Hep5, and Hep6 ( ). We also found through the KEGG enrichment analysis of Hep1/3/5/6 cells ( ) that the pathways involved in Hep1/3/5/6 cells mainly included glutathione metabolism, endocytosis, drug metabolism, and phagosomes, among which the common pathways were detoxification-related and endocytosis pathways. In addition, cell clusters of Hep5/6 cells also contained many fatty acid metabolism-related pathways. Of the 117 co-upregulated genes, 33.33% encoded ribosomal proteins, ( Down). Of the 39 co-downregulated genes, 71.79% were functionally annotated ( Down). To further understand the immune, detoxification, stress, and other related functions of the four clusters in the hepatopancreas, the key genes were selected for violin plot display before and after treatment. The detailed expression profiles of these genes in each cell population are clearly shown in ( ).

Figure 5

Enrichment analysis of differentially expressed genes in P. monodon hepatopancreatic cells under ammonia-stress. (A) Differential gene volcano map display of Hep1/3/5/6 clusters. (B) KEGG enrichment analysis of differential genes in Hep1/3/5/6 clusters. (C) The Venn diagram of coupregulated genes and its gene function in hepatopancreatic cells. (D) The Venn diagram of codownregulated genes and it gene function in hepatopancreatic cells. (E) Violin diagram showing functional genes in hepatopancreatic cells. (*p < 0.05 and **P < 0.01).

Discussion

In this study, high-throughput scRNA-seq based on the 10× Genomics platform was successfully performed on shrimp for the first time. The comprehensive scRNA-seq data provide information on the cell populations and their stress-activated states in both tissues. Importantly, our scRNA-seq could clearly distinguish 7 hepatopancreatic cell populations and 7 haemocyte populations and their developmental trajectories, which are less well understood and for which marker genes were not previously available.

Different haemocyte populations in the circulating haemolymph play different roles in crustacean immune responses, such as phagocytosis, encapsulation, nodule formation, and release of immune factors. Integrins, lysozyme C and ALF are important marker genes of Hem5. Integrins were considered as molecular marker candidates for shrimp phagocytosis (16). ALF is an important AMP in prawns, and multiple ALFs have been found to play important roles in the innate immune system of prawns (17). Lysozyme C is the only type of lysozyme known to be present in vertebrates, protozoa, and invertebrates, and is considered a molecular marker of some vertebrate macrophages (18). High expression of the tumour necrosis factor-α gene (TNF) is regarded as an important marker of mosquito “megacyte” cells (10). The TNF gene was also detected in the Hem5 cluster we identified, but it was not a significant marker gene for this cell group, which may also suggest shrimp haemocyte specificity. Therefore, we speculate that Hem5 may have functions similar to vertebrate macrophages; it is an important immune cell cluster for phagocytic activity and ALF production.

Compared with Hem5, Hem4 is a cell cluster that can secrete more types of antimicrobial peptides (AMPs). The high expression of AMPs (penaeidin, crustin, ALF, chelonianin-like), Peroxinectin and phenoloxidase-activating factors (PPAFs) are the main features of the Hem4 cell population. AMPs are natural defence molecules found in nearly all taxonomic kingdoms. Shrimp AMPs are mainly expressed in circulating haemocytes, the main site of the immune response (19). Peroxinectin has been confirmed to be a proPO system-related protein in P. monodon, and its adhesion activity does not require peroxidase activity (20). Peroxinectin is produced by the activation of the prophenoloxidase (proPO) system. Coincidentally, PPAFs are a group of serine proteases. They can convert prophenol oxidase (proPO) to the active form of phenol oxidase (PO) causing melanisation, which plays an important role in defence against invading pathogens (21). Coincidentally, PPAF, proPO, ALF and penaeidin are also considered to be important marker genes for haemocytes classification in shrimp Marsupenaeus japonicas (4). AMPs and melanisation are the two most powerful shrimp humoral responses. Intensive studies have demonstrated a close link between the proPO activation cascade and the AMP synthesis pathway via pattern recognition protein crosstalk in shrimp, which enables the innate immune system to execute an effective immune response (22, 23). Molecular processes such as peptidase activity and peptidase inhibitor activity are significantly enriched in Hem4, molecular processes related to cellular activities, such as chromosome organization and DNA replication are enriched in Hem6; cytoskeleton organization and apoptosis are enriched in Hem7. This finding suggests that Hem4 is the main effector cell population of humoral immunity characterized by antimicrobial peptide secretion and proPO system activation, while Hem6 and Hem7 are mainly involved in the activities of cells in the cellular responses, including apoptosis, phagocytosis, nodule formation, and encapsulation.

The hepatopancreas of decapod crustaceans is the central metabolic organ and a major target organ in response to environmental toxicants or pathogen stimuli. The cell population of Hep1 was identified based on high expression of multiple marker genes, including perlucin-like protein (PLP), peritrophin protein, lysozyme, trypsin and L-like procathepsin. PLP is a typical C-type lectin. It has recently been shown to bind to and agglutinate bacteria, affect bacterial phagocytosis and AMP expression, and play an important role in the innate immune response of shrimp Litopenaeus vannamei (24). Peritrophins are the major proteins that form the peritrophic membrane, which is an important component of the digestive and immune defence systems of shrimp and many insects (25). Studies in shrimp Metapenaeus ensis indicate that cathepsin L (CatL) is localized to hepatopancreatic digested B cells, whereas CatL mRNA is present in F cells but not in mature B cells. CatL is thought to be transcribed in one type of cell that rapidly evolves to a morphologically different cell with enzymatic functions (26). Interestingly, pseudotime analysis showed that Hep1 is the initiation point of development in the two cell development pathways in the hepatopancreas ( ). Therefore, we speculate that Hep1 cells are an important site for the synthesis of various digestive enzymes and lectins, and along with the formation of enzymes or proteins, they can develop into other types of cells to perform secretion, digestion and immune functions.

Hep3 and Hep5 in our scRNA-seq data are enriched in several genes encoding antimicrobial peptides (ALF, penaeidin), which are known as widely distributed immune effector molecules. Several genes that encode glutathione S-transferase (GST) of metabolic enzymes were significantly enriched in Hep2 cell population. As an important part of the body’s detoxification system, GST is responsible for catalysing the combination of GSH and xenobiotics, expelling glutathione conjugates from the body, and playing a role in detoxification and antioxidation (27). Moreover, in the pseudotime analysis, Hep2 was also significantly different from other clusters. In Drosophila single-cell sequencing studies, GST has been used as a molecular marker of cells to distinguish cell groups (9). The enrichment of pathways such as oxidative phosphorylation, TCA cycle, and glutathione metabolism further prove that Hep2 is a cell cluster involved in oxidative stress and detoxification metabolism. Overall, these findings lay a foundation for the detailed mapping of gene expression in shrimp cells and tissues, providing a comprehensive reference atlas of cellular diversity. The large number of marker genes identified in this study will facilitate future studies of crustacean cell diversity and function.

Cell types are dynamic in nature. Our scRNA-seq data provides a framework to distinguish cell types and their possible developmental routes and states. Haemocytes developed along Paths I and II, while the hepatopancreatic cells developed along Paths I, II, and III ( ). Monocle2 analysis showed that the expression levels of AMPs and GSH-Px showed different or even opposite trends at the developmental endpoints of haemocytes with different fates. In several invertebrates, AMPs are used as indicator genes for cell classification and cell state (4, 22). This finding may also suggest that AMPs and GSH-Px play important roles in cell differentiation and functional plasticity of shrimp.

In addition, hematopoietic cells have been described in some crustacean decapods, and the absence of proPO in hematopoietic cells is an important feature that distinguishes them from circulating blood cells (3, 28). The Hem2 cell population was detected expressing proliferation-related genes and cyclin-dependent kinase, but not proPO. Therefore, we propose Hem2 as a developmental initiation site with oligopotent potential. It is speculated that Hem2 can give rise to terminally-differentiated cell types and other possible states of activation, with Hem4 and Hem5 at the developmental end being the most mature immune cell types in haemocytes. The differentiation pathways of crustacean blood cells and hepatopancreatic cells are currently unclear, and our data provide new insights into this question.

The understanding of the immune response and detoxification process of ammonia-N stress will help to improve the stress resistance of shrimp, which is crucial for aquaculture production and environmental sustainability. A large number of genes and molecular pathways related to metabolism, detoxification, and immune defence are activated in different cell groups after ammonia-N stress ( , ). The gene expression of proPO in haemocytes (e.g., Hem3 and Hem4) increased rapidly, showing a positive feedback mechanism response to ammonia stress. However, the gene expression of proPO appears to be repressed in ammonia-stressed hepatopancreatic cells (e.g., Hep1, Hep3, Hep5). Invertebrates, including shrimp, have no acquired immunity. Melanisation, which is performed by activating the proPO system, plays an important role in the invertebrate innate immune system. Recent evidence has shown that blood cells with high proPO gene expression are specific to the melanisation responses, whereas cells with low proPO levels perform a broader range of biological tasks in mosquito (28). The differential expression of proPO between shrimp cell populations that we observed may represent distinct cell lineages and may also reflect functional specificity and diversity. The activation of proPO system in shrimp blood cells and the inactivity of hepatopancreatic cells may indicate that haemocytes are more involved in melanisation of the immune defence under ammonia-N stress. After ammonia-N treatment, genes such as insulin-degrading enzymes and steroid hormone proteins in haemocytes and GST in hepatopancreatic cells were significantly changed, and were mainly enriched in molecular pathways such as oxidative phosphorylation, glutathione metabolism, and phagosome. This finding revealed the response mechanism of antioxidation, detoxification, hormone metabolism, and phagocytosis after ammonia-N stress (7, 8, 29).

AMP, as an important defence line of innate immunity, also showed dramatic differences in expression after ammonia-N stress. The genes encoding a variety of AMPs showed a consistent expression trend in hepatopancreas and blood cells, and all of them were down-regulated in several cell clusters (e.g., Hem3, Hem4, Hep1, Hep3 and Hep6). Several studies have pointed out that excess ammonia-N suppresses the immune function of shrimp and increases the risk of microbial pathogen invasion or viral infection such as such as Vibrio alginolyticus, Lactococcus garveyi and white spot syndrome virus (30–32). Our results also support the above view, speculating that the inhibition of antimicrobial peptide expression and synthesis in blood cells and hepatopancreas under ammonia-N stress may be an important reason for the severe weakening of immune defence. It is worth mentioning that a large proportion (21.8% in blood and 33.3% in hepatopancreas) of ribosomal genes was significantly up-regulated in cells after ammonia-N stress. Recent work has shown that the synthesis of many defence proteins is required in stress response, and that ribosomal translation is crucial for maintaining timely synthesis of stress defence proteins (33). However, the effect of ribosomal genes on stress response has been largely overlooked in previous studies on shrimp. This lack of research also suggests that ribosome translation should be studied in parallel with traditional transcriptional regulation in the future to better understand shrimp stress defences.

Conclusion

By successfully applying scRNA-seq in shrimp, seven cell clusters were successfully identified in each of the two tissues based on specifically expressed marker genes. Our data reveal that Hem2 has oligopotent potential as the initiation site for haemocytes development. Hem4 is the main effector cell population of humoral immunity characterized by antimicrobial peptide secretion and activation of the proPO system. Hep1 cells are an important site for the synthesis of various digestive enzymes and lectins, and they can develop into other types of cells. A large number of DEGs are involved in oxidative stress, antioxidant system and proPO system, detoxification metabolism and immune defence. Several important genes including AMPs, proPO, and GST can be used as markers to identify cell populations, cell differentiation, and cell fate. Altogether, our scRNA-seq analysis documents the diversity of shrimp cell populations and provides a comprehensive resource of gene expression profiles for understanding the functions of cells under ammonia stress.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI SRA BioProject, accession no: PRJNA818766.

Ethics Statement

All experiments involving animals were approved by the Animal Care and Use Committee of South China Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences (No. SCSFRI96-253).

Author Contributions

SGJ and FZ conceived the project and supervised the work. YL performed the bioinformatics analysis and prepared the manuscript, tables and figures. ZM and JH conducted the experiment. LY, SJ, and QY collected the samples and performed sequencing. All authors read and approved the final manuscript.

Funding

This work was supported by the National Key R & D Program of China (2018YFD0900303), Research and Development Projects in Key Areas of Guangdong Province (2021B0202003), Central Public-interest Scientific Institution Basal Research Fund, CAFS (2020TD30), Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0605), China Agriculture Research System of MOF and MARA (CARS-48), Central Public Interest Scientific Institution Basal Research Fund, South China Sea Fisheries Research Institute, CAFS (2020ZD01, 2021SD13), Hainan Provincial Naturl Science Foundation of China (320QN359), Guangdong Basic and Applied Basic Research Foundation (2020A1515110200), and Guangzhou Science and Technology Planning Project (202102020208).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.