Abstract Simple Summary The Chinese soft-shelled turtle (Trionyx sinensis) is an important cultured reptile in East Asia. Hemorrhagic sepsis caused by Aeromonas hydrophila infection is the dominant disease in the aquaculture of Chinese soft-shelled turtles, while the molecular pathology is far from clear due to the lag of research on turtle immunology. It has been reported in mammals and fish that the dysfunction of immune responses to pathogen infections causes host tissue hemorrhagic sepsis. In this study, two groups of turtles with different susceptibility to A. hydrophila infection are found. A comparative transcriptome strategy is adopted to examine the gene expression profiles in liver and spleen for these two phenotypes of turtles post A. hydrophila infection, for the first time revealing the full picture of immune mechanisms against A. hydrophila, which provides new insight into the molecular pathology during A. hydrophila infection in T. sinensis. The findings will promote further investigations on pathogenic mechanisms of hemorrhagic sepsis caused by A. hydrophila infection in T. sinensis, and also will benefit their culture industry. Abstract Although hemorrhagic sepsis caused by Aeromonas hydrophila infection is the dominant disease in the aquaculture of Chinese soft-shelled turtle, information on its molecular pathology is seriously limited. In this study, ninety turtles intraperitoneally injected with A. hydrophila exhibited two different phenotypes based on the pathological symptoms, referred to as active and inactive turtles. Comparative transcriptomes of liver and spleen from these two groups at 6, 24, and 72 h post-injection (hpi) were further analyzed. The results showed that cytokine–cytokine receptor interaction, PRRs mediated signaling pathway, apoptosis, and phagocytosis enriched in active and inactive turtles were significantly different. Pro-inflammatory cytokines, the TLR signaling pathway, NLR signaling pathway, and RLR signaling pathway mediating cytokine expression, and apoptosis-related genes, were significantly up-regulated in inactive turtles at the early stage (6 hpi). The significant up-regulation of phagocytosis-related genes occurred at 24 hpi in inactive turtles and relatively lagged behind those in active turtles. The anti-inflammatory cytokine, IL10, was significantly up-regulated during the tested periods (6, 24, and 72 hpi) in active turtles. These findings offer valuable information for the understanding of molecular immunopathogenesis after A. hydrophila infection, and facilitate further investigations on strategies against hemorrhagic sepsis in Chinese soft-shelled turtle T. sinensis. Keywords: Chinese soft-shelled turtle, Aeromonas hydrophila, hemorrhagic sepsis, molecular immunopathogenesis 1. Introduction Chinese soft-shelled turtle (Trionyx sinensis) is an important reptile in East Asia and has been taken as a food resource for a long time in these areas, especially in China and Japan. In ancient Chinese medicine descriptions, consumption of turtles may bring about positive effects for human health, including strengthening immunity, anti-aging, and curing cardiocerebrovascular diseases [[34]1]. The excellent nutriment and medical values mean the turtles have become one of the most critical freshwater aquaculture reptiles, and the annual production is over 325,000 tons in China. However, infectious diseases caused by bacteria and viruses have resulted in severe losses to the aquaculture of turtles [[35]2,[36]3]. It is noteworthy that there are more than 15 kinds of diseases caused by Aeromonas hydrophila infection, such as red neck disease, septicemia, furunculosis, etc., accounting for about 60% of the total disease cases in turtles [[37]4,[38]5]. The histopathology of Chinese soft-shelled turtles infected by A. hydrophila has been previously described. In general, the pathogenic processes undergo the adhesion of bacteria adhesion factors, and the destruction of host liver, respiratory, and digestive organs by virulence factors such as exotoxin and extracellular enzymes, which finally leads to serious tissue hemorrhagic sepsis and the death of turtles [[39]6,[40]7]. Nevertheless, the molecular pathology of hemorrhagic sepsis caused by A. hydrophila infection is far from being elucidated in Chinese soft-shelled turtles. It has been reported in mammals and fish that the dysfunction of host immune responses to pathogen infections causes tissue hemorrhagic sepsis [[41]8,[42]9,[43]10]. The immune system generally serves as the security guard for hosts, which holds an immune network composed of many immune cells [[44]11]. The cells of the innate immune system consisting of neutrophils, monocytes, macrophages, dendritic cells, and natural killer cells recognize pathogens, produce cytokines, and engulf pathogens through phagocytosis, which are the first line of host defense against pathogens [[45]12]. These innate immune cells rely on pattern recognition receptors (PRRs) including toll-like receptors (TLRs), nucleotide-binding domain (NOD)-like receptors (NLRs), and retinoic acid-inducible gene I-like receptors (RLRs) to recognize various microbial invaders and produce cytokines that further activate the innate as well as adaptive immune cells [[46]13]. Under normal immune response, immune cells moderately produce and release cytokines, including interleukins (ILs), chemokines, and tumor necrosis factors (TNFs), which can contact more immune cells to participate in the battle between the host and foreign pathogens [[47]14]. When immune cells produce excessive cytokines, especially pro-inflammatory cytokines such as IL1, IL6, IL18, and TNFs (“cytokine storm”), the host immune system is overactivated and attacks self-tissues or cells, causing systemic inflammation, tissue hemorrhagic sepsis, and even death [[48]15,[49]16]. Although the research of turtle immunology lags behind that of mammals and fish, the unique evolutionary status as secondary aquatic reptiles has recently aroused wide concern on the distinct immune response mechanism against pathogens in turtles [[50]17,[51]18,[52]19,[53]20]. Fifteen candidate TLR family genes have been identified in T. sinensis [[54]21]. After A. hydrophila infection, TLR2 and TLR4 are significantly up-regulated in the spleen, indicating the immune response of TLR signaling pathway during bacterial infection in T. sinensis [[55]21]. Zhou et al. identified an IL8 homolog from T. sinensis and confirmed that IL8 mRNA expression shows significant up-regulation in various tissues, including liver, spleen, kidney, heart, intestine, and blood, after A. hydrophila infection [[56]22]. Zhang et al. also investigated the mRNA expression changes of pro-inflammatory cytokines such as IL1β, TNFα, IL6, IL8, and IL12 in T. sinensis during acute cold stress, revealing that acute cold stress increases the expression of pro-inflammatory cytokines in the spleen and intestine to withstand A. hydrophila infection [[57]23]. These above studies may help us understand the mRNA expression profiles and function of several immune molecules in T. sinensis. However, the evidence on the immune response mechanism in T. sinensis so far is too scattered and limited, which largely hinders the studying of pathogen–host immunity interaction and the molecular pathology during A. hydrophila infection in Chinese soft-shelled turtles. RNA-sequencing (RNA-Seq) using next-generation sequencing is one of the most useful methods to survey the character of a transcriptome because it offers the whole data on gene expression. To date, transcriptome profiling using next-generation sequencing technologies has provided new insights into pathogen–host immunity interaction in many aquaculture fish species, such as Nile tilapia (Oreochromis niloticus) [[58]24], Gold fish (Carassius auratus L.) [[59]25], Ussuri catfish (Pseudobagrus ussuriensis) [[60]26], Yellow catfish (Pelteobagrus fulvidraco) [[61]27], and Atlantic salmon (Salmo salar) [[62]28]. More and more studies have used this high-throughput sequencing approach to identify the expression differences of immune molecules between the resistant and susceptible fish to pathogen infection. For example, Moraleda et al. analyzed the comparative transcriptomes of resistant/susceptible salmons with different immune responses to Piscirickettsia salmonis infection and revealed that the gene networks involved in the apoptotic, cytoskeletal reorganization, bacterial invasion, and intracellular trafficking processes are tightly associated with disease resistance/susceptibility [[63]29]. We are originally inspired by the previous research results in mammals and fish that the dysfunction of host immune responses to pathogen infections causes tissue hemorrhagic sepsis. In the present study, therefore, a comparative transcriptome strategy is adopted to reveal the immune gene expression profiles in two phenotypes of Chinese soft-shelled turtles with different susceptibility to A. hydrophila infection at different time periods (6, 24, and 72 hpi), and to clarify which immune process abnormalities may be related to the occurrence of hemorrhagic sepsis. The results depict the full picture of immune mechanisms in response to A. hydrophila infection, suggesting that the dysfunction of cytokine–cytokine receptor interaction, PRRs mediated signaling pathways, phagocytosis, and apoptosis may cause hemorrhagic sepsis during A. hydrophila infection in Chinese soft-shelled turtles. This study, for the first time, reveals the host immunity–pathogen infection interaction and molecular immunopathogenesis during A. hydrophila infection by comparative transcriptomes, which may contribute to the development of novel management strategies for disease control and prevention in Chinese soft-shelled turtles. 2. Materials and Methods 2.1. Experimental Animals and Bacteria Strain Experimental turtles (T. sinensis) with an average weight of 16.45 ± 1.28 g were obtained from a breeding farm (111°97′ E, 28°90′ N) in Changde City, Hunan Province, China. All the turtles were acclimated in 50 L aquarium with aerated freshwater in a constant temperature laboratory room at 30 °C for two weeks before processing and were fed with a commercial diet (Kesheng Feed Stock Co., Ltd., Hangzhou, Zhejiang, China) twice a day. The animal experiments were according to the rules of the Animal Care and Use Committee of Hunan Agricultural University (Changsha, China; Approval Code: 201903297; Approval Date: 11 October 2019). The bacteria A. hydrophila, isolated from the spleen of clinically diseased turtles, was provided by Professor Zhipeng Gao from the College of Animal Science and Technology, Hunan Agricultural University, China. A. hydrophila was cultured in Luria Bertani (LB) medium at 30 °C with the shaking at 200 rpm. After 18 h culturing, the bacteria were harvested by centrifugation, and suspended with sterile phosphate buffer saline (PBS). The bacteria concentration was adjusted to 1.39 × 10^9 CFU/mL with PBS, which was previously proved as the medium lethal concentration for turtles during A. hydrophila infection. 2.2. Experimental Treatments, Pathological Observation, and Sampling For the A. hydrophila challenge, 90 turtles were intraperitoneally injected with bacteria suspension (100 μL per turtle, designated as the treatment group). After bacteria injection, the turtles in the treatment group could be divided into two subgroups based on their pathological symptoms and behavior activities. One subgroup of turtles showed weaker ability to feed and move, with significant hemorrhagic symptoms on the body surface and swelling and congestion on viscera after anatomy, which was designated as the inactive subgroup. The other subgroup of turtles exhibited no obvious pathological symptoms, which was designated as the active subgroup. Ten turtles, injected with 100 μL of PBS, were set up as the control group. To investigate the expressional profiles of immune-related genes during A. hydrophila infection in active and inactive turtles, important immune-related tissues including spleen and liver from the active subgroup turtles (N = 3) and inactive subgroup turtles (N = 3) were sampled post-injection of A. hydrophila at 6, 24, 72 h, respectively. Spleen and liver tissues from the control group turtles (N = 3) with PBS injection at 0 h were sampled as the control for the gene expression comparison with those from active subgroup turtles or from inactive subgroup turtles. Each sample was taken three times, and the sampled tissues were stored in liquid nitrogen before RNA extraction. 2.3. RNA Extraction, Library Preparation, and RNA Sequencing (RNA-Seq) Total RNA from each sampled tissue was isolated by using RNeasy mini kit (QIAGEN, Germantown, MD, USA) and treated with RNase-free DNase I (QIAGEN, Germantown, MD, USA) at 37 °C for 1 h to remove residual genomic DNA. RNA quality and concentration were determined using Agilent Bioanalyser 2100 (Agilent Technologies, Santa Clara, CA, USA) and NanoDrop-1000 (NanoDrop Technologies, Wilmington, DE, USA), respectively. Five micrograms of RNA were used to construct RNA-seq library according to the instruction of Illumina mRNA-Seq Prep Kit (Illumina, San Diego, CA, USA), and the libraries were sequenced by paired-end sequencing on the Illumina HiSeq 2500 sequencing platform (Illumina, San Diego, CA, USA). The quality of RNA-Seq raw reads were assessed with FastQC (version 0.10.1; [64]http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/ accessed on 11 March 2019), and were cleaned by removing adapter sequences, poly-N sequences, and low-quality sequences. The clean reads were then aligned to the published T. sinensis genome using HISAT2 [[65]30]. In total, 291.16 GB clean data were produced, and the average of 20.4~36.3 GB clean reads were obtained from samples. About 80% of the clean reads were mapped to the reference genome by using StringTie software (Center for Computational Biology, Johns Hopkins University, Baltimore, MD, USA) [[66]31], and the Q30 was >92.86%. The general information of RNA-Seq data is listed in [67]Table S1. 2.4. Identification of the Differentially Expressed Genes (DEGs) The relative transcript abundances in tissues (liver and spleen) from active and inactive subgroup turtles at different time periods of A. hydrophila infection compared to the control turtles were, respectively, estimated by using StringTie software with expectation maximization method, based on fragments per kilobase of exon per million fragments mapped (FPKM) [[68]31]. Differential expression analysis was performed by using the DESeq R package with default parameters [[69]32]. Benjamini and Hochberg’s approach was used to control the false discovery rate (FDR) by resulting P-values adjustment [[70]33]. Genes with an adjusted P-value (or q-value) < 0.05 found by DEGSeq were assigned as the differentially expressed. 2.5. KEGG Pathway Enrichment Analysis Significantly enriched signal transduction pathways represented by DEGs were determined using KEGG pathway enrichment analysis, compared with the whole genome background [[71]34]. The statistical enrichment of DEGs in KEGG pathways was tested using the software KOBAS, with a P-value < 0.05 [[72]34]. The significantly enriched KEGG pathways are listed in [73]Tables S2–S5. 2.6. Gene Expression Validation Using Quantitative PCR (qPCR) To validate RNA-seq data and gene expression profiles, nine DEGs were randomly selected to perform qPCR. Specific primers were designed using Primer 5 based on the coding sequences of identified genes from the turtle genome; the sequences of primers are listed in [74]Table S6. Primer specificity was ascertained using the following steps: PCR amplification, sequencing of PCR products, and BLAST analysis in the NCBI database. The RNA samples used for qPCR amplifications were the same as those used to construct the RNA-Seq library mentioned above. The qPCR was performed using the LightCycler^®480 Real-Time PCR System (Roche, Basel, Switzerland) with SYBR Green I Master. The reaction mixture (10 μL) comprised 2.5 μL cDNA, 0.5 μL (10 nM) forward primer and 0.5 μL (10 nM) reverse primer, 1.5 μL PCR grade water, and 5 μL Master Mix. Each reaction was performed in triplicate under the following conditions: 95 °C for 10 min, 40 cycles of 95 °C for 15 s, 55 °C for 15 s, and 72 °C for 30 s. The relative expression level of each gene was calculated according to the 2^−△△Ct method [[75]35] and normalized to the endogenous control genes of β-actin and GAPDH. 3. Results 3.1. Symptom Description of the Turtles Challenged with A. hydrophila Ninety Chinese soft-shelled turtles were injected with A. hydrophila (1.39 × 10^9 CFU/mL), and the pathological symptoms of turtles were observed. There were no obvious pathological symptoms on the body surface of turtles, and only six turtles showed behavior abnormality with slow-moving action at 6 hpi (post-injection of 6 h). At 24 hpi, the food intake and movement of turtles was overall reduced, and seven turtles died (accounting for about 8%), with the pathological symptom of abdominal congestion. At 72 hpi, 34 turtles (accounting for about 38%) died. The pathological symptoms on the body surface of diseased turtles were obvious, with white spots near the axillae, swelling, and congestion in chest and abdomen ([76]Figure 1A,B). After anatomizing the diseased turtles, pathological symptoms of turtle viscera were easily observed ([77]Figure 1C). The liver and spleen were swelling and congestive, with the tendency to decay, and the color of liver was yellow ([78]Figure 1C). The intestines were filled with no food debris, and the color was white ([79]Figure 1C). There were 13 turtles (accounting for about 14%) without any obvious pathological symptoms at 72 hpi ([80]Figure 1D). During the whole experiment within 72 h, 10 turtles in the control group showed no obvious pathological symptoms and were aggressive and active in feeding and moving (data not shown). Figure 1. [81]Figure 1 [82]Open in a new tab Pathological symptoms of turtles infected by A. hydrophila. The pathological symptoms on the body surface including (A) white spots near the axillae; (B) abdominal congestion in inactive subgroup turtles. (C) The pathological symptoms of viscera including liver, spleen, and intestines after anatomy of inactive turtles. (D) The active subgroup turtles showed no obvious pathological symptoms and were aggressive and active in feeding and moving. Black arrows indicate the location of pathological symptoms. According to pathological symptoms, the Chinese soft-shelled turtles challenged with A. hydrophila could be divided into two subgroups, including the active and inactive turtles. The liver and spleen were important immune organs in turtles, and both carried obviously pathological symptoms after A. hydrophila infection. Therefore, liver and spleen samples from active and inactive subgroup turtles were taken at 6, 24, and 72 hpi, respectively, and transcriptome sequencing was performed to reveal the differences of gene expression profiles between active and inactive turtles infected by A. hydrophila. Pearson’s correlation coefficients were first used to test for biologically repeated correlations between samples ([83]Figure 2A). The generated cluster dendrogram was used to observe the overall correlation of the transcriptomes from the AL group (liver in active turtles), the IL group (liver in inactive turtles), the AS group (spleen in active turtles), and the IS group (spleen in inactive turtles) at different time periods (6, 24, and 72 hpi) ([84]Figure 2A). Three biological replicates of liver and spleen samples from each time period and the transcriptome data both exhibited good correlation ([85]Figure 2A). The similarity test between the three biological replicates required the use of a principal component analysis (PCA) ([86]Figure 2B). Using the first principal component (PC1) and second principal component (PC2), a dimensionality reduction analysis was used to analyze the similarity between each replicate ([87]Figure 2B). [88]Figure 2B showed that biological replicates of samples overall exhibited good similarity. The generated box plot presented the dispersion degree of the gene expression level in a single liver or spleen sample, and intuitively revealed the whole gene expression level difference among all samples ([89]Figure 2C). The results showed that the gene expression level in the spleen was overall higher than that in the liver ([90]Figure 2C). Figure 2. [91]Figure 2 [92]Open in a new tab Transcriptional relationship between samples and the overall gene expression profiles. (A) Heatmap of correlation value (R square) of 37 libraries from liver or spleen samples. (B) Principal component analysis based on all of the expressed genes, showing 14 distinct groups of samples. (C) The dispersion degree of the gene expression level in a single liver or spleen sample. (D) The significantly up-regulated and down-regulated DEGs identified in livers or spleens from active and inactive turtles at 6, 24, and 72 hpi compared to the control. Totals of 4092 and 5793 DEGs were obtained in the liver and spleen transcriptomes, respectively ([93]Figure 2D). In the AL group, 321, 401, and 378 genes were significantly up-regulated, and 671, 874, and 151 genes were significantly down-regulated at 6, 24, and 72 hpi, respectively ([94]Figure 2D). In the IL group, the numbers of significantly up-regulated genes were 298, 138, and 79; the numbers of significantly down-regulated genes were 268, 381, and 132 at 6, 24, and 72 hpi, respectively ([95]Figure 2D). In the AS group, 219, 176, and 65 genes were significantly up-regulated, and 576, 437, and 206 genes were significantly down-regulated at 6, 24, and 72 hpi, respectively ([96]Figure 2D). In the IS group, the numbers of significantly up-regulated genes were 1193, 203, and 117; the numbers of significantly down-regulated genes were 1996, 401, and 204 at 6, 24, and 72 hpi, respectively ([97]Figure 2D). The results revealed that the number of DEGs in the spleen transcriptomes was more than that in the liver transcriptomes, and the molecular response peaked at 24 hpi in liver, while it peaked at 6 hpi in spleen ([98]Figure 2D). 3.2. Functional Classification of DEGs in Turtle Liver Transcriptomes by KEGG In order to investigate the different molecular response mechanisms against A. hydrophila infection in livers from active and inactive turtles, the functional classification of DEGs in AL and IL group transcriptomes at different time periods (6, 24, and 72 hpi) were analyzed by KEGG enrichment analysis, and the results are summarized as follows. 3.2.1. Sequential Changes of KEGG Enrichment in AL Group Turtles In AL group, the up-regulated DEGs were mainly enriched in immune-related pathways including “cytokine–cytokine receptor interaction”, phagocytosis-associated processes including “phagosome”, “protein processing in endoplasmic reticulum”, and “protein export”, and pathogen infection-related pathways, while the down-regulated DEGs were intensively involved in 17 metabolism pathways and three cell adhesion-related processes at 6 hpi ([99]Figure 3). Figure 3. [100]Figure 3 [101]Open in a new tab KEGG enrichment analysis in AL group turtles at different time periods. The top 20 KEGG pathways are presented here in the form of scatterplots to show the up-regulated and down-regulated DEGs enriched in livers from active subgroup turtles at 6, 12, and 72 hpi. The enrichment factor is the ratio between the DEG number and the number of all genes in a certain gene enrichment term. The sizes of the dots on these plots denote the number of DEGs, while colors correspond to the q value range. At 24 hpi, the majority of up-regulated DEGs were also annotated into immune-related pathways including “cytokine–cytokine receptor interaction”, “phagosome”, “protein processing in endoplasmic reticulum”, “proteasome”, and “protein export” ([102]Figure 3). In addition, cytokine expression-mediating pathways including “toll-like receptor signaling pathway” and “NOD-like receptor signaling pathway” were significantly up-regulated ([103]Figure 3). Similar to the response at 6 hpi, the down-regulated DEGs were mainly involved in a series of metabolism and cell adhesion pathways ([104]Figure 3). At 72 hpi, the up-regulated DEGs could be functionally classified into “toll-like receptor signaling pathway” and “RIG-I-like receptor signaling pathway”, pathogen infection-related pathways, “apoptosis”, and several metabolism-related pathways ([105]Figure 3), while the down-regulated DEGs mainly participated in important metabolism-related pathways such as “insulin signaling pathway” and “adipocytokine signaling pathway” ([106]Figure 3). 3.2.2. Sequential Changes of KEGG Enrichment in IL Group Turtles Unlike the KEGG enrichment in AL group at 6 hpi, the up-regulated DEGs were mainly enriched in cytokine expression-mediating pathways including “toll-like receptor signaling pathway”, “NOD-like receptor signaling pathway”, and “RIG-I-like receptor signaling pathway” and “apoptosis”, besides “cytokine–cytokine receptor interaction”, pathogen infection-related pathways; while phagocytosis-associated processes were not listed in the top 20 up-regulated KEGG pathways in IL group at 6 hpi ([107]Figure 4). The down-regulated DEGs were functionally annotated into cell adhesion-related pathways such as “ECM–receptor interaction” and “focal adhesion”, and energetic metabolism pathways ([108]Figure 4). Figure 4. [109]Figure 4 [110]Open in a new tab KEGG enrichment analysis in IL group turtles at different time periods. The top 20 KEGG pathways are presented here in the form of scatterplots to show the up-regulated and down-regulated DEGs enriched in livers from inactive subgroup turtles at 6, 12, and 72 hpi. The enrichment factor is the ratio between the DEG number and the number of all genes in a certain gene enrichment term. The sizes of the dots on these plots denote the number of DEGs, while colors correspond to the q value range. At 24 hpi, the term of “cytokine–cytokine receptor interaction” was not listed in the top 20 up-regulated KEGG pathways, while the up-regulated DEGs were mainly involved in phagocytosis-associated processes such as “phagosome”, “protein processing in endoplasmic reticulum”, “proteasome”, and “protein export” ([111]Figure 4). Additionally, “toll-like receptor signaling pathway” was significantly up-regulated ([112]Figure 4). For the down-regulated DEGs, most of them were associated with hormone synthesis and amino acid metabolism ([113]Figure 4). At 72 hpi, immune-related terms including “toll-like receptor signaling pathway” and “cytosolic DNA-sensing pathway” were listed in the top 20 up-regulated KEGG pathways ([114]Figure 4). The down-regulated DEGs were mainly enriched in energetic metabolism pathways such as “insulin signaling pathway” and “adipocytokine signaling pathway” ([115]Figure 4). 3.2.3. Expression Difference Analysis of Cytokine, Phagocytosis, and Apoptosis-Related Genes between AL and IL Group Turtles The KEGG pathways related to immune processes including cytokine–cytokine receptor interaction, phagocytosis, and apoptosis enriched in AL and IL group turtles challenged with A. hydrophila were quite different ([116]Table 1, [117]Table 2 and [118]Table 3). Therefore, the fold changes of differentially expressed cytokine, phagocytosis, and apoptosis-related genes were further analyzed. Table 1. Fold changes of differentially expressed cytokine and cytokine receptor, and cytokine expression mediating pathway genes in the AL group and the IL group compared to the control. Categories/ Gene Name Description Log[2] (Fold Changes) AL Group IL Group 6 Hpi 24 Hpi 72 Hpi 6 Hpi 24 Hpi 72 Hpi Interleukins and interleukin receptors IL1β interleukin 1 beta 2.18 4.82 IL8 interleukin 8 3.55 IL10 interleukin 10 6.08 4.41 2.31 IL1R1 interleukin 1 receptor type I 1.66 IL1R2 interleukin 1 receptor type II 5.88 3.16 7.06 IL5RA interleukin 5 receptor alpha −2.57 −2.40 −3.01 IL18R1 interleukin 18 receptor 1 −5.30 Chemokines and chemokine receptors CCL5 C-C motif chemokine 5 4.71 CCL20 C-C motif chemokine 20 9.34 7.63 12.18 CX3CL1 C-X3-C motif chemokine 1 −4.91 3.28 CCR5 C-C chemokine receptor type 5 3.02 TNF family members and TNF receptors TNFSF10 tumor necrosis factor ligand superfamily member 10 −2.56 −2.08 TNFSF15 tumor necrosis factor ligand superfamily member 15 −2.40 SF6B tumor necrosis factor receptor superfamily member 6B −2.57 SF12A tumor necrosis factor receptor superfamily member 12A 4.49 4.60 Toll-like receptor (TLR) signaling pathway TLR2 toll-like receptor 2 4.56 2.80 TLR5 toll-like receptor 5 3.61 1.51 2.89 TRIF toll-like receptor adapter molecule 1 −3.15 −2.87 −1.99 −3.40 −3.42 TAB1 TAK1-binding protein 1 1.77 PI3K phosphoinositide-3-kinase 2.48 2.96 MAP2K1 mitogen-activated protein kinase kinase 1 2.13 MAP2K6 mitogen-activated protein kinase kinase 6 −2.04 -2.63 −2.59 AP-1 proto-oncogene protein c-fos 1.43 STAT1 signal transducer and activator of transcription 1 1.49 RIG-I-like receptor (RLR) signaling pathway RIG-I retinoic acid inducible gene I 2.38 LGP2 laboratory of genetics and physiology 2 2.21 MDA5 melanoma differentiation-associated gene 5 2.43 TRAF3 TNF receptor-associated factor 3 5.07 3.54 IRF7 interferon regulatory factor 7 −1.83 2.86 DDX3X ATP-dependent RNA helicase DDX3X −1.86 −1.69 −2.08 −2.02 NOD-like receptor (NLR) signaling pathway ASC apoptosis-associated speck-like protein containing a CARD 1.99 RIPK2 receptor-interacting serine/threonine-protein kinase 2 2.47 cIAP baculoviral IAP repeat-containing protein 2/3 2.50 TNFAIP3 tumor necrosis factor alpha-induced protein 3 1.56 3.75 [119]Open in a new tab Note: This table and the following tables only show the DEGs with the values of log[2] (fold change) with FDR < 0.05. Table 2. Fold changes of differentially expressed phagocytosis-related genes in the AL group and the IL group compared to the control. Categories/ Gene Name Description Log[2] (Fold Changes) AL Group IL Group 6 Hpi 24 Hpi 72 Hpi 6 Hpi 24 Hpi 72 Hpi Internalization and formation of the phagosomes TLR2 toll-like receptor 2 4.56 2.80 MR mannose receptor 1.87 1.07 −2.18 iC3b the fragment of complement component 3 4.38 2.79 Collectin C-type lectin 4.81 3.69 1.47 −2.63 F-actin actin beta/gamma 1 5.57 5.31 4.39 5.09 5.68 4.64 Early phagosome Rab5 ras-related protein Rab-5B 2.28 vATPase V-type H^+-transporting ATPase 2.33 2.05 3.85 −2.08 CALR calreticulin 4.05 Mature phagosome TUBA tubulin alpha 2.26 TUBB tubulin beta 4.63 3.25 1.42 2.72 vATPase V-type H^+-transporting ATPase 2.33 2.05 3.85 −2.08 Phagolysosome sec61 protein transport protein SEC61 subunit beta 2.72 3.26 2.95 vATPase V-type H^+-transporting ATPase 2.33 2.05 3.85 −2.08 Activation of NADPH oxidase p40phox neutrophil cytosolic factor 4 3.50 p47phox neutrophil cytosolic factor 1 2.79 2.42 gp91 NADPH oxidase 1 1.59 Antigen presentation MHC II MHC class II antigen −2.21 sec22 vesicle transport protein SEC22 4.21 3.54 [120]Open in a new tab Table 3. Fold changes of up-regulated apoptosis-related genes in the AL group and the IL group compared to the control. Gene Name Description Log[2] (Fold Changes) AL Group IL Group 6 Hpi 24 Hpi 72 Hpi 6 Hpi 24 Hpi 72 Hpi p53 tumor protein p53 1.94 IP3R inositol 1,4,5-triphosphate receptor type 3 1.35 1.56 1.92 Perforin perforin 1 3.12 2.57 PI3K phosphoinositide-3-kinase 2.48 2.96 MERK2 mitogen-activated protein kinase kinase 2 2.13 PERK protein kinase RNA (PKR)-like ER kinase 1.85 2.30 Cathepsin cathepsin B 5.11 4.27 NOXA phorbol-12-myristate-13-acetate-induced protein 1 1.62 AP1 proto-oncogene protein c-fos 1.43 GZMB granzyme B 3.87 IL3R cytokine receptor common subunit beta 3.43 A1 hematopoietic Bcl-2-related protein A1 4.59 [121]Open in a new tab The results showed that the sharpest response of cytokines occurred at 24 hpi in the AL group, with the up-regulation of two cytokine and two cytokine receptor genes ([122]Table 1). Among three cytokine expression-mediating pathways, only the response of toll-like receptor signaling pathway was relatively intense, with five up-regulated genes, also at 24 hpi in the AL group ([123]Table 1), while in IL group, the sharpest response of cytokines occurred at 6 hpi, with the up-regulation of five cytokine and two cytokine receptor genes, and few cytokines were differentially expressed at 24 or 72 hpi ([124]Table 1). Moreover, the up-regulation of differentially expressed cytokine genes at 6 hpi in IL group was overall higher than those in AL group at any time periods ([125]Table 1). Nevertheless, IL10, an anti-inflammatory cytokine, was not differentially expressed at any time periods in the IL group, while it was up-regulated at all the time periods in the AL group ([126]Table 1). In addition, the up-regulations of both NOD-like receptor signaling pathway and RIG-I-like receptor signaling pathway genes were intense in the IL group at 6 hpi ([127]Table 1). Notably, the sharpest response of phagocytosis-related gene occurred at 6 hpi in the AL group, with the up-regulation of 10 DEGs, and lasted until 24 hpi ([128]Table 2). While in the IL group, at 24 hpi, the phagocytosis activity seemed to just start, with five up-regulated phagocytosis-related DEGs ([129]Table 2). In addition, the most intense response of apoptosis-related genes occurred at 6 hpi in the IL group, with the up-regulation of six DEGs, while in the AL group, the apoptosis seemed to just start at 24 hpi with the up-regulation of five DEGs ([130]Table 3). 3.3. Functional Classification of DEGs in Turtle Spleen Transcriptomes by KEGG The functional classification of DEGs in AS and IS transcriptomes at different time periods were also analyzed by KEGG enrichment to further reveal the different immune response mechanisms against A. hydrophila infection in spleens between active and inactive turtles. The results are summarized as follows. 3.3.1. Sequential Changes of KEGG Enrichment in AS Group Turtles In the AS group, the up-regulated DEGs were functionally associated with immune processes including “cytokine–cytokine receptor interaction”, “chemokine signaling pathway”, “phagosome”, “apoptosis”, “leukocyte transendothelial migration”, “toll-like receptor signaling pathway”, “MAPK signaling pathway”, and pathogen infection-related pathways at 6 hpi ([131]Figure 5), while the down-regulated DEGs were involved in cell adhesion and metabolism pathways at 6 hpi ([132]Figure 5). Figure 5. [133]Figure 5 [134]Open in a new tab KEGG enrichment analysis in AS group turtles at different time periods. The top 20 KEGG pathways are presented here in the form of scatterplots to show the up-regulated and down-regulated DEGs enriched in spleens from active subgroup turtles at 6, 12, and 72 hpi. The enrichment factor is the ratio between the DEG number and the number of all genes in a certain gene enrichment term. The sizes of the dots on these plots denote the number of DEGs, while colors correspond to the q value range. At 24 hpi, the up-regulated DEGs were mainly annotated into phagocytosis-related pathways such as “phagosome”, “protein processing in endoplasmic reticulum”, “synaptic vesicle cycle”, and “lysosome”. In addition, “toll-like receptor signaling pathway” was significantly up-regulated ([135]Figure 5). Unlike the response at 6 hpi, the down-regulated DEGs at 24 hpi could be enriched in pathways including “cytokine–cytokine receptor interaction” and “bacterial invasion of epithelial cells”, besides a series of metabolism and cell adhesion pathways ([136]Figure 5). At 72 hpi, the up-regulated DEGs were functionally classified into apoptosis-related pathways, including “p53 signaling pathway” and “cell cycle”, as well as immune defense processes such as “phagosome”, “chemokine signaling pathway”, and “leukocyte transendothelial migration” ([137]Figure 5), while the down-regulated DEGs mainly participated in important metabolism and cell adhesion pathways ([138]Figure 5). 3.3.2. Sequential Changes of KEGG Enrichment in IS Group Turtles The strongest immune response of spleen occurred at 6 hpi in the IS group, with 1193 up-regulated and 1996 down-regulated DEGs. Most up-regulated DEGs were enriched in immune-related pathways including “cytokine–cytokine receptor interaction”, “toll-like receptor signaling pathway”, “NOD-like receptor signaling pathway”, and “RIG-I-like receptor signaling pathway”, as well as apoptosis-associated processes ([139]Figure 6), while the down-regulated DEGs were functionally annotated into cell adhesion-related pathways such as “ECM–receptor interaction” and “focal adhesion”, and metabolism pathways ([140]Figure 6). Figure 6. [141]Figure 6 [142]Open in a new tab KEGG enrichment analysis in IS group turtles at different time periods. The top 20 KEGG pathways are presented here in the form of scatterplots to show the up-regulated and down-regulated DEGs enriched in spleens from inactive subgroup turtles at 6, 12, and 72 hpi. The enrichment factor is the ratio between the DEG number and the number of all genes in a certain gene enrichment term. The sizes of the dots on these plots denote the number of DEGs, while colors correspond to the q value range. At 24 hpi, the up-regulated DEGs could be related to phagocytosis such as “phagosome” and “lysosome” ([143]Figure 6). Additionally, “toll-like receptor signaling pathway” was listed in the top 20 up-regulated KEGG pathways ([144]Figure 6). For the down-regulated DEGs, most of them were associated with cell adhesion, hormone synthesis, and amino acid metabolism ([145]Figure 6). Similar to the response at 72 hpi in the AS group, in the IS group, the up-regulated DEGs mainly participated in apoptosis-related pathways including “p53 signaling pathway” and “cell cycle”, as well as “phagosome” and “chemokine signaling pathway” at 72 hpi ([146]Figure 6), while the down-regulated DEGs were enriched in a series of hormone synthesis, amino acid metabolism, and cell adhesion pathways ([147]Figure 6). 3.3.3. Expression Difference Analysis of Cytokine, Phagocytosis, and Apoptosis-Related Genes between AS and IS Group Turtles The fold changes of differentially expressed cytokine, phagocytosis, and apoptosis-related genes were also analyzed in AS and IS group turtles ([148]Table 4, [149]Table 5 and [150]Table 6). The results showed that in both AS and IS groups, the sharpest response of cytokines occurred at 6 hpi, while the number of differentially expressed cytokine and cytokine receptor genes in the IS group were overwhelmingly more than that in the AL group (30 up-regulated genes in the IS group vs. seven up-regulated genes in the AS group) ([151]Table 4). Noteworthily, the overall expression of cytokines dramatically decreased at 24 hpi and 72 hpi in the IS group, while gradually decreased in the AS group, and the anti-inflammatory cytokine IL10 was up-regulated at all the time periods in the AS group ([152]Table 4). Consistent with the response of cytokines, the expression of toll-like receptor signaling pathway genes were also relatively intense at 6 hpi, but only with three up-regulated genes in the AS group ([153]Table 4), while in the IS group at 6 hpi, the expressions of three cytokine expression-mediating pathway genes, including toll-like receptor signaling pathway, NOD-like receptor signaling pathway, and RIG-I-like receptor signaling pathway, were overall up-regulated, with 18 up-regulated genes ([154]Table 4). Table 4. Fold changes of differentially expressed cytokine and cytokine receptor, and cytokine expression mediating pathway genes in the AS group and the IS group compared to the control. Categories/ Gene Name Description Log[2] (Fold Changes) AS Group IS Group 6 Hpi 24 Hpi 72 Hpi 6 Hpi 24 Hpi 72 Hpi Interleukins and interleukin receptors IL1β interleukin 1 beta 2.46 IL6 interleukin 6 5.20 IL7 interleukin 7 −2.79 IL8 interleukin 8 3.27 2.43 2.41 IL10 interleukin 10 6.82 5.69 1.21 7.30 IL1R2 interleukin 1 receptor type II 4.13 3.81 IL1RAP interleukin 1 receptor accessory protein 2.24 IL3RB cytokine receptor common subunit beta 3.56 IL4R interleukin 4 receptor 1.18 IL5RA interleukin 5 receptor alpha −1.95 −3.81 −2.85 IL5RB interleukin 5 receptor alpha 3.56 IL8RB interleukin 8 receptor 3.07 2.32 IL12RB1 interleukin 12 receptor beta-1 2.60 IL12RB2 interleukin 12 receptor beta-1 3.20 IL15RA interleukin 15 receptor alpha 3.68 IL21R interleukin 21 receptor 1.98 IL22RA2 interleukin 22 receptor alpha 2 5.21 4.77 4.39 Chemokines and chemokine receptors CCL20 C-C motif chemokine 20 6.33 6.72 2.56 CXCL10 C-X-C motif chemokine 10 5.36 CXCL11 C-X-C motif chemokine 11 4.32 CXCL12 C-X-C motif chemokine 12 −2.00 CXCL13 C-X-C motif chemokine 13 3.47 2.41 CXCL14 C-X-C motif chemokine 14 1.81 2.07 CX3CL1 C-X3-C motif chemokine 1 −3.70 3.90 CCR5 C-C chemokine receptor type 5 1.26 CXCR4 C-X-C chemokine receptor type 4 −1.75 XCR1 XC chemokine receptor 1 −4.81 TNF family members and TNF receptors TNFSF8 tumor necrosis factor ligand superfamily member 8 3.06 TNFSF10 tumor necrosis factor ligand superfamily member 10 −2.86 −4.23 −3.17 TNFSF12 tumor necrosis factor ligand superfamily member 12 −1.99 TNFSF15 tumor necrosis factor ligand superfamily member 15 3.09 TNFSF18 tumor necrosis factor ligand superfamily member 18 4.56 EDA ectodysplasin-A −2.93 SF6B tumor necrosis factor receptor superfamily member 6B 7.29 SF9 tumor necrosis factor receptor superfamily member 9 3.26 SF12A tumor necrosis factor receptor superfamily member 12A 6.32 7.32 7.37 SF13B tumor necrosis factor receptor superfamily member 13B 2.33 3.08 SF19L tumor necrosis factor receptor superfamily member 19-like 5.98 FAS tumor necrosis factor receptor superfamily member 6 2.34 NGFR tumor necrosis factor receptor superfamily member 16 −3.24 EDAR ectodysplasin-A receptor −4.72 Toll-like receptor (TLR) signaling pathway TLR4 toll-like receptor 4 2.30 1.71 TLR5 toll-like receptor 5 2.77 1.82 3.00 MyD88 myeloid differentiation factor 88 1.08 TRIF toll-like receptor adapter molecule 1 1.96 Rac Ras-related C3 botulinum toxin substrate 1 2.37 2.52 1.97 1.84 2.51 2.04 PI3K phosphoinositide-3-kinase −4.06 AKT RAC serine/threonine-protein kinase 1.26 MAP2K1 mitogen-activated protein kinase kinase 1 2.88 2.82 2.50 2.85 2.94 2.57 MAP2K6 mitogen-activated protein kinase kinase 6 −2.37 AP-1 proto-oncogene protein c-fos 1.60 1.33 STAT1 signal transducer and activator of transcription 1 2.03 IKBKE inhibitor of nuclear factor kappa-B kinase subunit epsilon 1.92 RIG-I-like receptor (RLR) signaling pathway RIG-I retinoic acid inducible gene I 2.32 LGP2 laboratory of genetics and physiology 2 1.90 MDA5 melanoma differentiation-associated gene 5 1.54 TRAF2 TNF receptor-associated factor 2 2.19 MITA Mediator of IFN regulatory transcription factor 3 activation 1.29 NOD-like receptor (NLR) signaling pathway NLRP12 NACHT, LRR and PYD domains-containing protein 12 −2.65 −2.05 −3.11 ASC apoptosis-associated speck-like protein containing a CARD 2.27 RIPK2 receptor-interacting serine/threonine-protein kinase 2 1.46 CARD8 caspase recruitment domain-containing protein 8 1.49 [155]Open in a new tab Table 5. Fold changes of differentially expressed phagocytosis-related genes in the AS group and the IS group compared to the control. Categories/ Gene Name Description Log[2] (Fold Changes) AS Group IS Group 6 Hpi 24 Hpi 72 Hpi 6 Hpi 24 Hpi 72 Hpi Internalization and formation of the phagosome TLR4 toll-like receptor 4 2.30 1.71 MR mannose receptor −2.71 −1.43 −3.02 CR1 complement receptor 1 1.69 2.15 αvβ5 integrin αvβ5 −2.45 iC3b the fragment of complement component 3 −3.36 Collectin C-type lectin −4.00 −1.66 −4.85 F-actin actin beta/gamma 1 3.60 −3.11 −3.66 −3.71 −2.46 Early phagosome Rab5 ras-related protein Rab-5B 2.47 1.84 1.82 vATPase V-type H^+-transporting ATPase 2.21 −1.01 1.59 CALR calreticulin 2.10 Mature phagosome TUBA tubulin alpha 1.14 1.58 1.63 TUBB tubulin beta 3.63 1.60 −1.21 2.28 2.34 vATPase V-type H^+-transporting ATPase 2.21 −1.01 1.59 Dynein dynein heavy chain 2 −1.82 Phagolysosome sec61 protein transport protein SEC61 subunit beta 1.36 vATPase V-type H^+-transporting ATPase 2.21 −1.01 1.59 NOS nitric-oxide synthase −3.07 TAP ATP-binding cassette subfamily B member 3 2.17 Activation of NADPH oxidase gp91 NADPH oxidase 1 1.62 1.57 1.40 p40phox neutrophil cytosolic factor 4 1.27 p47phox neutrophil cytosolic factor 1 2.75 2.32 2.53 2.04 1.52 P67phox neutrophil cytosolic factor 2 1.40 Rac Ras-related C3 botulinum toxin substrate 1 2.37 2.52 1.97 1.84 2.51 2.04 Antigen presentation MHC II MHC class II antigen 3.14 −1.87 −2.54 [156]Open in a new tab Table 6. Fold changes of up-regulated apoptosis-related genes in the AS group and the IS group compared to the control. Gene Name Description Log[2] (Fold Changes) AS Group IS Group 6 Hpi 24 Hpi 72 Hpi 6 Hpi 24 Hpi 72 Hpi p53 tumor protein p53 1.66 IP3R inositol 1,4,5-triphosphate receptor type 3 1.17 Perforin perforin 1 3.24 5.58 MERK2 mitogen-activated protein kinase kinase 2 2.88 2.82 2.50 2.85 2.95 2.57 Cathepsin cathepsin B 6.19 3.24 1.57 2.81 AP1 proto-oncogene protein c-fos 1.60 1.33 GZMB granzyme B 1.22 1.57 5.59 5.34 IL-3R cytokine receptor common subunit beta 2.56 TRAF12 TNF receptor-associated factor 2 2.91 2.19 Fas tumor necrosis factor receptor superfamily member 6 2.34 TrkA neurotrophic tyrosine kinase receptor type 1 6.49 NIK mitogen-activated protein kinase kinase kinase 14 1.24 FLIP CASP8 and FADD-like apoptosis regulator 1.26 eiF2α translation initiation factor 2 subunit 1 1.25 Calpain calpain-1 1.49 ARTS septin 4 2.45 AIF apoptosis-inducing factor 1 1.40 Gadd45 growth arrest and DNA-damage-inducible protein 2.32 ASK1 mitogen-activated protein kinase kinase kinase 5 5.27 CytC cytochrome c 1.90 1.80 [157]Open in a new tab The intense response of phagocytosis started at 6 hpi in the AS group, with the up-regulation of seven DEGs, and lasted until 24 hpi (six up-regulated DEGs) ([158]Table 5), while in the IS group at 6 hpi, the phagocytosis activity seemed to be inhibited, with 10 down-regulated DEGs, and just started at 24 hpi, with 11 up-regulated phagocytosis-related DEGs ([159]Table 5). Nevertheless, in both the AS and IS groups, the sharpest response of apoptosis-related gene expression occurred at 6 hpi, while the number of up-regulated apoptosis-related genes in the IS group (18 DEGs) was significantly more than that in the AL group (four DEGs) ([160]Table 6). 3.4. Validation of DEGs by qPCR A total of nine DEGs were randomly selected to perform qPCR to validate RNA-Seq data and gene expression profiles ([161]Figure 7). PCR products with expected sizes were successfully amplified with all the nine specific primer pairs, indicating their availabilities for DEG validation (data not shown). The different amplification efficiencies of the nine DEGs between the experimental and control groups were transformed by log[2] (fold change) to compare with the results of RNA-Seq. The results showed that the expression patterns of these genes determined by qPCR were similar to those acquired through RNA-Seq ([162]Figure 7), which confirmed the reliability of the RNA-Seq data. Therefore, the immune-related genes isolated in this study could be useful references for future studies on the molecular mechanisms of