Immune Response and Transcriptome Analysis of the Head Kidney to Different Concentrations of Aeromonas veronii in Common Carp (Cyprinus carpio)

Abstract

The common carp (Cyprinus carpio), a major economic freshwater fish, is suffering from a variety of bacterial infectious diseases because of its high-density, factory and intensive farming patterns. Aeromonas veronii is the causative agent of high mortality in common carp, causing severe economic losses in aquaculture. However, the regulatory mechanisms involved in the response of common carp to this bacterial infection remain poorly understood. In this study, we compared mortality, blood serum LZM (Lysozyme) and IgM (Immunoglobulin M) levels and transcriptome patterns of head kidney tissues after infection with different concentrations of Aeromonas veronii. We observed that mortality increased progressively with an increasing pathogen concentration. The concentrations of blood serum LZM and IgM significantly increased after infection. A total of 13 and 925 differentially expressed genes (DEGs) were identified after infection with low (T4) and high (T9) concentrations of bacterial suspension, respectively. KEGG and GO analyses of the DEGs highlighted multiple immune-related signaling pathways. Weighted gene co-expression network analysis (WGCNA) revealed that 136 and 83 hub genes were related to blood serum LZM and IgM, respectively. Finally, the gene expression in the head kidney was validated via RT–qPCR to be consistent with the transcriptome. These results provide insights into the mechanisms of the immune response to infection with different concentrations of Aeromonas veronii and offer useful information for further studies on immune defense mechanisms in common carp.

1. Introduction

The common carp (Cyprinus carpio) is a prominent freshwater species that is widely distributed throughout the world, thriving in a vast array of environmental conditions due to its remarkable adaptability [1,2,3]. The common carp has a long domestication and cultivation history in China and not only plays an important role in the traditional Chinese culture but also is a delicacy of people’s tables, especially in northern China [4,5]. External environmental stresses pose a threat to the healthy farming of common carp by weakening their immune system and increasing their sensitivity to various pathogens [6,7]. In addition, high-density artificial farming and deterioration of water quality exacerbated the rapid spread of diseases, resulting in large−scale economic losses. It is an urgent issue that needs to be addressed to reduce the risk of disease and improve the overall amount of healthy common carp farming.

Aeromonas veronii, Gram−negative parthenogenetic anaerobic bacilli, is a type of Aeromonas that is a critical pathogen in the aquaculture industry and is especially widespread in freshwater [8,9]. It is an opportunistic pathogen in fish that causes bacterial infectious diseases, such as septicemia, hemorrhage, ulcers and ascites, and high mortality in diseased fish, resulting in considerable economic losses in aquaculture and posing a serious risk to fisheries worldwide [10,11,12]. Notably, this bacterium can produce zoonotic diseases through aquatic products, and its harm as an aquatic zoonotic agent has been demonstrated by an increasing number of cases [13,14,15]. Many studies have focused on the mechanisms of resistance to Aeromonas veronii infection in fish. The pathogenesis of this bacterium depends mainly on its virulence genes and virulence factors. The strain SJ4 has been reported to cause cellular enlargement, marked hemorrhage and inflammatory responses in diseased Siniperca chuatsi, and SJ4 not only carried virulence genes such as act, fim, flgM, ompA, lip, hly, aer and eprCAL, but produced tyrosinase, dnase, gelatinase and hemolysin [16]. Changes in the expression levels of host-associated immune genes are key mechanisms involved in the response to the infection. Aeromonas veronii could cause mass mortalities in Odontobutis potamophila, and the immune-related gene expression levels of MHC II, Myd88, TLR and SOD were significantly increased in the liver, gill, spleen and head kidney at different infection times [17]. A total of 114 immune−related DEGs were captured in the spleen transcriptome profiles of African catfish (Clarias gariepinus) challenged with this bacterium and were significantly enriched in 38 pathways related to immunity or disease, including the NF-kappa B, TNF, NLR, TLR and RLR pathways [18]. In addition, key genes involved in resistance to its infection in certain fish species have been identified, such as galectin−9, which was involved in the immune response against infection in koi carp (Cyprinus carpio) [19]; the CD4-1 gene, which activated NF−κB signaling in response to infection [20]; and the mRNA expression levels of CaCD3γ/δ, which significantly changed in the spleen, head kidney, intestine and gill after infection in Qihe crucian carp (Carassius auratus) [21].

To explore the molecular mechanisms underlying the defense against Aeromonas veronii infection, we compared the transcriptome profile differences between head kidney tissues from control carp and infection groups with different concentrations of bacterial suspension and identified immune-related genes and signaling pathways activated or suppressed during this pathogen infection. Understanding the mechanisms of the immune response to Aeromonas veronii in common carp is crucial for controlling and preventing infections in aquaculture development strategies. Ultimately, this research facilitates the development of scientifically informed approaches to increase the disease resistance of freshwater aquaculture species.

2. Results

2.1. Comparison of Morbidity and Mortality at Six Infection Concentrations

Within 14 days of Aeromonas veronii infection, we counted morbid individuals, dead individuals and the time of death in the six groups with different infection concentrations. Except for the T4 group, the other five infection groups included morbid individuals with obvious symptoms, including skin ulceration, hemorrhages on the body surface, and intestinal dropsy (Figure 1a). In the T5 and T6 groups, only two morbid individuals were observed. Nevertheless, a total of 22 morbid individuals and three dead individuals were observed in the T7 group, and the morbidity and mortality were 73.33% and 10.00%, respectively. The highest morbidity of 96.67% was recorded in the T8 and T9 groups, and the mortality reached 80.00% in the T9 group. Overall, morbidity and mortality increased with increasing concentrations of Aeromonas veronii infection, especially in the high-concentration T9 group, compared with those in the low-concentration T4 group (Table 1 and Figure 1b). Notably, three acutely dead individuals in the T9 group died at 1 day post infection due to high concentrations of Aeromonas veronii (Supplementary Table S1).

2.2. Blood Serum LZM and IgM Responses to Different Infection Concentrations

To assess the immune response to Aeromonas veronii infection, the concentrations of blood serum LZM and IgM in common carp were compared across four representative groups with different concentrations of infection. We selected the samples in the low-concentration T4 group (n = 10) and the high-concentration T7 (n = 6), T8 (n = 10) and T9 (n = 6) groups. In the noninfected group, the blood serum LZM concentration was greater than the IgM concentration in the NC group, at 1040.95 ± 211.19 pg/mL and 12.85 ± 2.91 pg/mL, respectively (Supplementary Table S2). On the 14th day after infection, the concentrations of LZM and IgM in all infection groups were significantly greater than those in the NC group (p ≤ 0.05) (Figure 2). Among the infection groups, the blood serum LZM concentration increased with increasing infection concentrations, reaching a maximum level in the T8 group (1880.95 ± 269.87 pg/mL), but it decreased in the T9 group (1667.84 ± 238.59 pg/mL). The blood serum IgM concentration was highest in the T4 group (18.83 ± 2.59 pg/mL) and decreased with an increasing concentration, indicating that LZM and IgM exhibited different patterns of immune response when common carp were challenged with different concentrations of Aeromonas veronii.

2.3. Identification and Analysis of DEGs Under Low- and High-Concentration Infection

Since the expression patterns of key genes in head kidney tissues reflect the level of the immune response, we compared differences in the expression profiles of head kidney tissues in response to different concentrations of Aeromonas veronii infection. We obtained at least 42,341,064 clean reads from each sample. The average Q20 and Q30 values of the clean reads were 98.67% and 96.44%, respectively, indicating that the sequencing quality of all the samples was satisfactory (Supplementary Table S3).

Compared with those in the NC group, 13 DEGs were detected in the low-concentration T4 group, including five upregulated and eight downregulated genes (Supplementary Table S4), whereas 925 DEGs were detected in the high-concentration T9 group, including 478 upregulated and 447 downregulated genes (Supplementary Table S5) (Figure 3). The number of DEGs was greater in the high-concentration T9 group than in the low-concentration T4 group, suggesting that a higher concentration of Aeromonas veronii elicits a stronger immune response at the gene expression level. In particular, a total of four shared DEGs were detected in two pairs: T4 vs. NC and T9 vs. NC. In addition, nine and 921 unique DEGs remained in the T4 vs. NC and T9 vs. NC pairs, respectively.

2.4. Functional and Pathway Enrichment Analyses of the Identified DEGs

KEGG and GO pathway enrichment analyses were performed to provide a deeper understanding of the functions of the DEGs. A total of 13 DEGs and 925 DEGs obtained from T4 vs. NC and T9 vs. NC groups, respectively, were annotated with KEGG terms. No significant signaling pathways were enriched in the T4 vs. NC group. However, in the T9 vs. NC group, upregulated genes were significantly enriched in 15 KEGG terms, including the Toll−like receptor signaling pathway, the RIG−I−like receptor signaling pathway, the NOD−like receptor signaling pathway, Cytokine−cytokine receptor interaction, the C−type lectin receptor signaling pathway and other important immune signaling pathways (Figure 4b and Supplementary Table S6). Twelve significant KEGG terms were obtained from the downregulated genes, with metabolic pathways being the top significant. (corrected p value = 4.10 × 10⁻¹³) (Figure 4b and Supplementary Table S7). GO enrichment of the DEGs among the upregulated genes in the T9 vs. NC group revealed that 32 out of the top 41 enriched belonged to the biological process category (Figure 5a and Supplementary Table S8). Moreover, among the downregulated genes, 37 of the 64 GO terms belonged to the biological process category (Figure 5b and Supplementary Table S9). The positive regulation of response to stimulus, regulation of the immune system process, immune response, and immune system process were enriched in this category. These pathways play important roles in resistance to Aeromonas veronii infection in the head kidney tissues. Apparently, more immune signaling pathways were activated in head kidney tissues of common carp when they were exposed to high concentrations of Aeromonas veronii.

2.5. WGCNA Revealed the Modules Relevant to Serum LZM and IgM

WGCNA was performed to investigate the associations between gene modules and particular immune parameters in the serum. After filtering, a total of 15 out of 17 samples and 19,954 out of 43,531 genes were used to construct the co−expression network based on their expressional levels. A scale-free network was subsequently constructed on the basis of the soft threshold, and a clustering dendrogram was generated from gene clustering on topological overlap matrix (TOM) dissimilarity (Figure 6a,b). The correlations between module feature vectors and serum LZM and IgM were analyzed, and no modules were associated with LZM, but two modules were highly correlated with serum IgM: MEgreen (p = 0.04, correlation coefficient = 0.54) and MEpurple (p = 0.04, correlation coefficient = 0.53) (Figure 6c). Therefore, as modules of interest, the two modules were further studied below.

2.6. Functional and Pathway Enrichment Analysis of the Significant Modules and Screening of Hub Genes

A total of 2915 and 514 genes were co-expressed in the green and purple modules, respectively (Figure 6a,b). To elucidate the biological significance of these genes, KEGG pathway enrichment revealed that 24 and 15 terms reached a significant level (BH method) (p ≤ 0.05) in the two modules. Coincidentally, the cell cycle, DNA replication, proteasome, oocyte meiosis, RNA polymerase, RNA transport and Fanconi anemia pathways were coenriched in two modules (Figure 7c,d) (Supplementary Tables S10 and S11). We subsequently identified 136 hub genes from the 2915 genes in the green module (Supplementary Table S12) and 83 hub genes from the 514 genes in the purple module (Supplementary Table S13), according to gene significance (GS ≥ 0.2) and module membership (MM ≥ 0.8).

2.7. Validation of RNA-Seq Results

To further verify the expression patterns of the genes identified in the RNA−Seq expression analysis and hub genes, we randomly selected five genes for RT−qPCR analysis, namely, col6a3 (collagen, type VI, alpha 3), smad6b (SMAD family member 6b), dsg2.1 (desmoglein 2, tandem duplicate 1), LOC109069404 (collagen alpha−2(IV) chain−like) and zfp36l2 (zinc finger protein 36, C3H type−like 2), and four hub genes, LOC109053336 (C−C motif chemokine 8−like), LOC109048270 (SLAM family member 5-like), LOC109053842 (proteinase−activated receptor 1−like) and LOC109096249 (collagen alpha−1(IV) chain−like). The same expression trends were found for the col6a3, smad6b, dsg2.1, LOC109069404 and zfp36l2 genes in the two pairs. Those genes were downregulated in the treated groups. Only two hub genes, LOC109053336 and LOC109048270, were upregulated. In summary, the RT−qPCR results were consistent with the transcriptome sequencing results, demonstrating the accuracy of the transcriptome sequencing data (Figure 8a,b).

3. Discussion

The common carp is a farmed fish species worldwide and an important food source for people, indicating its outstanding economic importance in freshwater aquaculture. Despite their strong environmental adaptability, common carp are still exposed to a variety of pathogens, restricting the development of healthy farming. In particular, the spread of Aeromonas veronii is a great threat to the economy in freshwater aquaculture. However, the potential molecular mechanism of immune responses against this pathogen infection remains poorly understood. Some studies have reported the immune response to Aeromonas veronii infection in various fish species [22,23,24,25], but these studies are limited in their ability to explain the mechanisms of this pathogen infection in common carp. Our study provides insights into the immune defense mechanisms of common carp in response to infection with different concentrations of Aeromonas veronii.

Aeromonas veronii can cause high mortality in carp, resulting in obvious symptoms, such as skin ulceration, hemorrhages on the body surface and abdominal and intestinal dropsy in diseased fish, which was similar with previous reports in other freshwater fish [22,24]. In this study, we detected more obvious symptoms and higher mortality under challenge with high concentrations of bacterial suspension. In addition, the diseased fish exhibit a typical inflammatory response, as evidenced by a significant increase in the blood serum LZM and IgM levels. Changes in these two indicators at different infection concentrations have rarely been reported. Notably, LZM, a hydrolytic enzyme present in the body surface mucus, intestinal mucus, serum and macrophages of many fish species, plays an important role in the fish innate immune system and participates in resistance to pathogens such as bacteria, fungi and viruses. IgM is one of the few immunoglobulins found in fish to date and has powerful bactericidal and immunomodulatory effects. Among the infection groups, the trend of increased levels of these two immune factors in the blood serum of diseased fish was different, indicating that there was a time lag between different concentrations of bacteria inducing increases in LZM and IgM. For example, in the T4 group, the level of IgM was higher than that in the other groups at 14 days post-infection, and it is possible that a higher bacterial load triggered the immune response earlier, which coincided with the time of death of the fish in each group.

The head kidney is one of the most important lymphoid tissues in teleosts and plays a crucial role in the immune system against microbial pathogens [26,27]. Changes in gene expression in the head kidney constitute an intrinsic mechanism of response to pathogenic bacterial infection. We compared the DEGs of head kidney tissues subjected to low and high concentrations of Aeromonas veronii infection via transcriptome analysis. Compared with those in the control group, greater concentrations of pathogen infection caused more genes to be differentially expressed in head kidney tissues. Pathway enrichment analysis of the DEGs revealed that multiple immune signaling pathways were activated in the host defense system against pathogenic bacterial infection. KEGG and GO analyses revealed the enrichment of the Toll−like receptor signaling pathway, the RIG-I-like receptor signaling pathway, the NOD−like receptor signaling pathway, Cytokine−cytokine receptor interaction, the C−type lectin receptor signaling pathway, the positive regulation of response to stimulus, regulation of immune system process, immune response and immune system process, which was consistent with the results of Aeromonas veronii infection in other fish. For example, the Toll-like receptor signaling pathway and Cytokine−cytokine receptor interaction were enriched after Aeromonas veronii infection in channel catfish [28]. Transcriptome profiles of the spleen of African catfish (Clarias gariepinus) challenged with Aeromonas veronii revealed that NF−kappa B, TNF, NLR, TLR and RLR pathways related to immunity were significantly enriched [18]. In addition, similar immune signaling pathways were activated in Mastacembelus armatus during bacterial infection [29].

Changes in the expression of immune-related genes were key factors in regulating these immune signaling pathways. Among the verified genes, smad6b is a type of SMAD family member 6 and is involved in the regulation of inflammation. SMAD family member 6, a critical mediator of the TGF-β–BMP pathway, negatively regulates interleukin 1-receptor−Toll−like receptor signals through direct interaction with the adaptor Pellino-1 [30]. dsg2.1 encodes desmoglein 2, tandem duplicate 1, which is located mainly in the plasma membrane and is involved in homophilic cell adhesion via plasma membrane adhesion molecules [31]. The relationship between desmoglein 2 and inflammation has been reported, for example, intracellular desmoglein 2 cleavage sensitizes epithelial cells to apoptosis in response to proinflammatory cytokines [32]. zfp36l2 is a zinc finger protein, which is zinc finger protein 36, C3H type−like 2. Zinc finger proteins play important roles in gene regulation. Several zinc finger proteins are crucial for immune regulation because they target mRNAs for degradation and the modulation of signaling pathways [33]. There are no direct reports that the above genes are associated with Aeromonas veronii infection in common carp, but our data revealed that the above genes were downregulated in the head kidney after Aeromonas veronii infection, suggesting a consistent pattern of the immune response of these genes in the head kidney in response to infection with different concentrations of Aeromonas veronii. WGCNA was performed to further identify immune-related genes. Among the hub genes, LOC109053336, encoding a C−C motif chemokine 8-like protein, is an inflammatory mediator capable of recruiting immune cells at the site of infection. As a kind of kinase−interacting protein, akip1 plays a key role in immunological processes. Gao et al. reported that Akip1 acts as a molecular determinant of cAMP−dependent protein kinase in the NF−kappaB signaling pathway, which plays a crucial role in the pathogenesis of many NF−kappaB−related diseases [34]. LOC109048270 encodes the signaling lymphocytic activation molecule (SLAM) family member 5−like protein, and SLAM family receptors mediate important regulatory signals in immune cells by binding to members of the SLAM−associated protein (SAP) family of adaptors [35,36].

The host immune response to invading pathogenic bacteria is a systematic and complex biological process. Our study revealed the mechanism underlying the varying degrees of response of common carp to different infectious concentrations of Aeromonas veronii, including mortality, blood serum immune factors, and differentially expressed genes in head kidney tissues. Although we found no direct evidence that DEGs or hub genes regulated the immune response in common carp, the biological functions of these genes are closely related to the process of immunity, and further research on these genes is needed.

4. Materials and Methods

4.1. Sampling and Ethics Statement

This study complied with the regulations of the Animal Care and Use Committee of the Chinese Academy of Fishery Sciences (ACUC-CAFS). A panmictic population was cultured at the Fangshan experimental base of the Chinese Academy of Fishery Sciences (Beijing, China) and fed commercial feed (TongWei, Chengdu, China). A total of 210 half-aged juvenile fish with no skin damage were randomly selected for the experiment (mean body weight: 133.29 ± 26.76 g, mean body length: 17.80 ± 1.21 cm). Prior to the start of the experiment, all the fish were temporarily reared in recirculating water for 7 days at 20 °C.

4.2. Bacterial Strains and Infection

Aeromonas veronii was used as the infectious strain, which was recorded at China General Microbiological Culture Collection Center (CGMCC 7231). Sterile Luria–Bertani (LB) liquid and LB agarose gel mediums were prepared before bacterial rejuvenation. After the glycerol bacteria preserved at −80 °C were thawed, the tip of a sterilized pipette gun was dipped into the bacterial preservation suspension and then streaked on an LB agar medium and incubated at 28 °C overnight. A single colony was then picked, streaked on new LB plates, and incubated overnight at 28 °C. Then, a single colony was picked and inoculated into an LB liquid medium and incubated overnight at 28 °C with shaking. The subsequently obtained bacterial suspension was measured with an OD = 0.876. Based on the linear correlation between the OD value and concentration (correlation coefficient of 5 × 10⁸), the concentration of the bacterial suspension was calculated to be 4.38 × 10⁹ cfu/mL. The bacterial suspension was diluted 10 times with sterile saline medium. A total of 210 fish were randomly divided into seven groups, with 30 individuals per group. The control group (NC) consisted of 30 fish. The other 180 fish were infected with 4.38 × 10⁴ (T4), 4.38 × 10⁵ (T5), 4.38 × 10⁶ (T6), 4.38 × 10⁷ (T7), 4.38 × 10⁸ (T8), 4.38 × 10⁹ (T9) (cfu/mL) Aeromonas veronii, respectively. For each test group, 100 μL of bacterial suspension was injected into the back muscle of each fish. Each sample in the NC group was injected with the same volume of a sterile saline medium. Morbid fish visually exhibited signs of skin ulceration and hemorrhages on the body surface. The experiment lasted for 14 days, and the number of dead individuals in each group during the experiment was recorded. Head kidney tissues were taken from all surviving individuals after the experiment. Before the tissues were collected, all the fish were euthanized with chlorbutanol (0.05%).

4.3. Detection of Serum LZM and IgM

On the 14th day post infection, the concentrations of blood serum LZM and IgM in 42 surviving individuals from the NC (n = 10), T4 (n = 10), T7 (n = 6), T8 (n = 10) and T9 (n = 6) groups were measured via enzyme-linked immunosorbent assay (ELISA) (Bioswamp, Wuhan, China). In brief, standards and samples in duplicate were added to the appropriate Microelisa strip plate wells precoated with an antibody specific to LZM or IgM and combined with the specific antibody. Then, a horseradish peroxidase (HRP)-conjugated antibody specific for LZM or IgM was added to each MicroELISA strip plate well and incubated. After the free components were washed away, TMB substrate suspension was added to each well. Only those wells containing LZM or IgM and HRP-conjugated antibodies appeared blue and then turned yellow after the addition of the stop solution. The optical density (OD) was measured spectrophotometrically at a wavelength of 450 nm. The OD value was proportional to the concentration of LZM or IgM, and the concentrations of LZM or IgM in the samples were calculated by comparing the OD of the samples to the standard curve.

4.4. RNA Extraction

To compare the gene expression levels of head kidney tissues in the immune response to Aeromonas veronii, we selected the lowest (T4) and highest (T9) concentration groups as representative infected groups for comparative analysis with the NC group to identify differentially expressed genes (DEGs). A total of 17 head kidney tissue samples from the NC (n = 7), T4 (n = 7) and T9 (n = 3) groups were used in the present study. The total RNA of each sample was extracted via a TRIzol^® reagent (Invitrogen, Carlsbad, CA, USA). The concentration of RNA was measured with a Nanodrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA). RNA integrity was assessed with the RNA Nano 6000 Assay Kit of the Agilent 2100 bioanalyzer (Agilent Technologies, CA, USA). The RNA purity was determined via agarose gel electrophoresis.

4.5. RNA-Seq Library Construction and Sequencing

For each sample, an mRNA sequencing library with a size of approximately 250–300 bp was generated via an NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, MA, USA). The library was sequenced on an Illumina platform in 150 bp paired-end mode.

4.6. Quality Control and Annotation

We used Trimmomatic to remove the read pairs with adapters, reads containing N and low-quality reads (reads with Qphred ≤ 5 bases accounting for more than 50% of the entire length of the read) [37]. The clean reads were aligned to the reference genome of common carp (GCF_018340385.1) [38] with HISAT2 [39]. StringTie was used to estimate the gene expression levels and generate new transcripts [40]. The new transcripts were annotated by alignment to the Pfam [41] and KEGG [42] databases. The gene expression level was measured as fragments per kilobase of transcript per million fragments mapped (FPKM) [43].

4.7. DEGs and Pathway Enrichment Analysis

We identified the differentially expressed genes (DEGs) related to the immune response to Aeromonas veronii in T4 and T9 groups, respectively. DEseq2 was used to detect differentially expression genes (DEGs) [44]. The p values were corrected via the Benjamini–Hochberg (BH) method [45], and genes with the parameters of corrected p value ≤ 0.05 and |log2(fold change)| ≥ 1 were assigned as DEGs. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was applied to identify the functions of the DEGs via KOBAS [46]. We used KOBAS to compare common carp proteins with zebrafish proteins and obtained Gene Ontology (GO) annotations of all the genes. The GO terms of DEGs were obtained with TBtools [47]. A term with a corrected p value ≤ 0.05 was considered to be significant.

4.8. Gene Co-Expression Network Analysis

The transcriptomes were used to construct a co−expression network via the R package “WGCNA” [48] to identify gene modules closely related to the serum LZM and IgM. The soft−thresholding parameter β was set to 6 on the basis of the criterion of approximate scale−free topology R² ≥ 0.8. The weighted adjacency matrix was transformed into a topological overlap matrix (TOM) to generate networks according to the network connectivity of each gene [49,50,51]. To group genes with similar expression profiles into the same gene modules, we performed average association hierarchical clustering on the basis of TOM-based similarity (1−TOM) with a minimum size cut−off of 100 and a deep split value of 2 for the resulting dendrogram. Dendrites representing highly relevant genes were grouped into a module and labelled with a specific color, and unassigned background genes that did not belong to any module were indicated by gray sections. Module eigengenes (MEs) were referred to as “hub genes” and dependent on degree of module membership (MM) and gene significance (GS), which were the cut-off criteria (|MM| ≥ 0.8 and |GS| ≥ 0.2) [52].

4.9. Validating the Gene Expression via Real-Time Quantitative PCR (RT–qPCR)

The genes identified in the RNA−Seq expression analysis and hub genes were randomly selected and validated via RT–qPCR. Total RNA from the head kidney tissues of another nine individuals, including the NC (n = 3), T4 (n = 3) and T9 (n = 3) groups, was extracted with a total RNA rapid extraction kit (BIOFIT, China). cDNA was synthesized with a FastKing one−step degenomic cDNA first strand synthesis premix reagent (TIANGEN, Beijing, China). The specific primers (Table 2) were designed via Primer-BLAST tools at the National Center for Biotechnology Information (NCBI). β−actin was used as an internal reference gene [53]. RT–qPCR was conducted on a Life Real 9600 thermal cycler using PowerUpTM SYBRTMGreen Master Mix (Applied Biosystems, CA, USA). The thermal cycling conditions were 95 °C for 3 min and 40 cycles of two steps (95 °C for 3 s and 60 °C for 34 s). The relative expression levels of genes between the NC group and the bacterial infection T4 and T9 groups were calculated via the 2^−ΔΔCT method [54].

4.10. Statistical Analysis

The significant differences between any two groups were determined with one-way analysis of variance (one−way ANOVA). All the results were visualized via R−Project (version 4.2.1) or GraphPad Prism software version 8.0. The data are expressed as the means ± standard deviations (SDs) unless otherwise stated.

5. Conclusions

In conclusion, our study revealed that high concentrations of Aeromonas veronii infection can cause high mortality and significantly increase the levels of blood serum LZM and IgM in common carp. The transcriptome profiles revealed that the DEGs were involved mainly in immune-related pathways, including the Toll−like receptor signaling pathway, the RIG−I−like receptor signaling pathway, the NOD−like receptor signaling pathway, Cytokine–cytokine receptor interaction, the C−type lectin receptor signaling pathway, the positive regulation of response to stimulus, regulation of immune system process, immune response and immune system processes. Furthermore, we identified 136 and 83 hub genes from the immune-related green and purple modules, respectively. These results revealed the mechanism of the immune response to Aeromonas veronii infection in common carp, providing a theoretical foundation for research on disease resistance in fish.

Author Contributions

J.Z.: Writing−original draft, Writing−review and editing, visualization, software, data curation. N.D.: Validation, data curation. Y.Q.: Validation, data curation. N.J.; Resources, conceptualization. W.X.: Resources. T.L.: Resources. Z.M.: Resources. Y.C.: Resources, investigation. Y.Z.: Supervision, resources, project administration, methodology. J.L.: Supervision, resources, project administration, methodology. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

This study was conducted in accordance with the Animal Care and Use Committee of the Chinese Academy of Fishery Sciences (ACUC−CAFS) under the established recommendations for the care and use of animals in scientific research.

Informed Consent Statement

Informed consent was obtained from all subjects involved in this study.

Data Availability Statement

The transcriptome data reported in this research are available at the Sequence Read Archive (SRA) under accession number PRJNA1128011.

Acknowledgments

We would like to thank Lin Luo for her technical support.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of this study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Symptoms of diseased fish and survival curves of common carp in the six groups challenged with different concentrations of Aeromonas veronii. (a) Symptoms of diseased fish: the left image shows the skin ulceration phenotype, and the right image shows the abdominal anatomy; (b) Survival curves of common carp in the six groups challenged with different concentrations of Aeromonas veronii. The four lines of NC, T4, T5 and T6 coincided.

Figure 2. Comparison of the serum LZM and IgM concentrations in common carp challenged with different infection concentrations. (a) Concentration differences of LZM in different groups; (b) Concentration differences of IgM in different groups. The five groups of data with different letters indicate significant differences in the figure (p ≤ 0.05).

Figure 3. DEGs in the low-concentration and high-concentration groups compared with those in the NC group. (a) Statistics of DEGs in two pairs: T4 vs. NC and T9 vs. NC. (b) Intersection analysis of differentially expressed genes between T4 vs. NC and T9 vs. NC. (c) Heatmap of DEGs in T4 vs. NC: the horizontal coordinate represents the sample number. (d) Heatmap of DEGs in the T9 vs. NC: the horizontal coordinate represents the sample number.

Figure 4. The enriched KEGG terms obtained from KEGG pathway enrichment analysis of the DEGs. (a) KEGG terms from enrichment analysis of upregulated genes in T9 vs. NC. (b) KEGG terms from enrichment analysis of downregulated genes in T9 vs. NC.

Figure 5. The enriched GO terms obtained from the DEGs. (a) GO terms from enrichment analysis of upregulated genes in T9 vs. NC. (b) GO terms from enrichment analysis of downregulated genes in T9 vs. NC.

Figure 6. Gene co-expression network analysis of serum LZM and IgM. (a) Scale-free topology model and mean connectivity. (b) Clustering dendrogram of genes, with dissimilarity based on the topological overlap, together with assigned module colors. The clustered branches represent different modules, and each line represents one gene. (c) Module associations with serum LZM and IgM levels, two values in each box; the upper values represent the absolute values of the correlation coefficients, and the lower values represent the p values.

Figure 7. KEGG analysis of genes in the green and purple modules. (a) Scatterplot of gene significance and module membership in the green module. The upper right corner represents the hub genes (GS ≥ 0.2 and MM ≥ 0.8). (b) Scatterplot of gene significance and module membership in the purple module. The upper right corner represents the hub genes (GS ≥ 0.2 and MM ≥ 0.8). (c) Pathway enrichment analysis of all genes in the green module. (d) Pathway enrichment analysis of all genes in the purple module.

Figure 8. Validation of RNA-Seq results via RT−qPCR. (a) Comparison of the results of transcriptome sequencing and RT−qPCR between the T4 and NC groups. (b) Comparison of the results of transcriptome sequencing and RT−qPCR between the T9 and NC groups.

Table 1

Morbidity and mortality statistics for the six groups with different pathogen concentrations.

Group	Concentration (cfu/mL)	Total Number	Morbid Individuals	Dead Individuals	Morbidity	Mortality
NC	0	30	0	0	0.00%	0.00%
T4	4.38 × 104	30	0	0	0.00%	0.00%
T5	4.38 × 105	30	2	0	6.67%	0.00%
T6	4.38 × 106	30	2	0	6.67%	0.00%
T7	4.38 × 107	30	22	3	73.33%	10.00%
T8	4.38 × 108	30	29	9	96.67%	30.00%
T9	4.38 × 109	30	29	24	96.67%	80.00%

Table 2

Primer sequences for RT–qPCR.

Gene Name	Forward Primer (5′-3′)	Reverse Primer (5′-3′)
β-actin	GATGATGAAATTGCCGCACTG	ACCAACCATGACACCCTGATGT
col6a3	GCCATCCGTGAGTTTATTCGG	CTTCCACCTTTTGCTGTAAGGC
smad6b	GACTTTCCCCCAGCGAAGAA	CGCCACGTTACACCAATGAC
dsg2.1	CTTCAAAGACGGCAACAGGC	TCCCCATAGTTCCTTCTGGGA
zfp36l2	GTTTGCTCACGGCTATCACG	GAGCCTCGGCTGAGTGTTAC
LOC109069404	GGGGACCAGAGGGACCTAAT	CCTGGTGGGCCTTTGAAGAA
LOC109096249	TCAAGATGGCTCTCCAGGTTTA	CAAATCTAGCAGGCTCTCCCT
LOC109053842	TGTTTCCCATCGCCTCTCTG	TGTCATTTTTGCGTTGCACAT
LOC109053336	GAGCTGCTGTAAGGAGGTGT	CGGTCTTTAATCCAGCGCAC
LOC109048270	GTGTTGGAGGGAGGTGTTCTC	GTCTAACTGCACTCTGCCTGT

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms252212070/s1.

References

1. Koehn, J.D. Carp (Cyprinus carpio) as a powerful invader in Australian waterways. Freshw. Biol.; 2004; 49, pp. 882-894. [DOI: https://dx.doi.org/10.1111/j.1365-2427.2004.01232.x]

2. Schilling, H.T.; Butler, G.L.; Cheshire, K.J.M.; Gilligan, D.M.; Stocks, J.R.; Thiem, J.D.; Crook, D.A. Contribution of invasive carp (Cyprinus carpio) to fish biomass in rivers of the Murray–Darling Basin, Australia. Biol. Invasions; 2024; 26, pp. 2955-2971. [DOI: https://dx.doi.org/10.1007/s10530-024-03362-x]

3. Fanson, B.G.; Hale, R.; Thiem, J.D.; Lyon, J.P.; Koehn, J.D.; Bennett, A.F.; Stuart, I. Assessing impacts of a notorious invader (common carp Cyprinus carpio) on Australia’s aquatic ecosystems: Coupling abundance-impact relationships with a spatial biomass model. Biol. Conserv.; 2024; 290, 110420. [DOI: https://dx.doi.org/10.1016/j.biocon.2023.110420]

4. Nakajima, T.; Hudson, M.J.; Uchiyama, J.; Makibayashi, K.; Zhang, J. Common carp aquaculture in Neolithic China dates back 8000 years. Nat. Ecol. Evol.; 2019; 3, pp. 1415-1418. [DOI: https://dx.doi.org/10.1038/s41559-019-0974-3]

5. Yao, G.; Li, W. Mandarin Fish Culture: Status and Development Prospects: Success Stories and Modern Trends. Aquaculture in China; Gui, J.F.; Tang, Q.S.; Li, Z.J.; Liu, J.S.; Silva, S.S. John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2018; pp. 256-269.

6. Zivna, D.; Sehonova, P.; Plhalova, L.; Marsalek, P.; Blahova, J.; Prokes, M.; Divisova, L.; Stancova, V.; Dobsikova, R.; Tichy, F. et al. Effect of salicylic acid on early life stages of common carp (Cyprinus carpio). Environ. Toxicol. Pharmacol.; 2015; 40, pp. 319-325. [DOI: https://dx.doi.org/10.1016/j.etap.2015.06.018]

7. Zhou, M.; Qiang, J.; Gan, J.; Xu, X.; Li, X.; Zhang, S.; Xu, B.; Dong, Z. Quercetin attenuates environmental Avermectin-induced ROS accumulation and alleviates gill damage in carp through activation of the Nrf2 pathway. Comp. Biochem. Physiol. C Toxicol. Pharmacol.; 2023; 274, 109744. [DOI: https://dx.doi.org/10.1016/j.cbpc.2023.109744]

8. Parte, A.C. LPSN–list of prokaryotic names with standing in nomenclature. Nucleic Acids Res.; 2014; 42, pp. D613-D616. [DOI: https://dx.doi.org/10.1093/nar/gkt1111]

9. Abbott, S.L.; Cheung, W.K.; Janda, J.M. The genus Aeromonas: Biochemical characteristics, atypical reactions, and phenotypic identification schemes. J. Clin. Microbiol.; 2003; 41, pp. 2348-2357. [DOI: https://dx.doi.org/10.1128/JCM.41.6.2348-2357.2003]

10. Rahman, M.; Colque-Navarro, P.; Kühn, I.; Huys, G.; Swings, J.; Möllby, R. Identification and characterization of pathogenic Aeromonas veronii biovar sobria associated with epizootic ulcerative syndrome in fish in Bangladesh. Appl. Environ. Microbiol.; 2002; 68, pp. 650-655. [DOI: https://dx.doi.org/10.1128/AEM.68.2.650-655.2002]

11. Vega-Sánchez, V.; Acosta-Dibarrat, J.; Vega-Castillo, F.; Castro-Escarpulli, G.; Aguilera-Arreola, M.G.; Soriano-Vargas, E. Phenotypical characteristics, genetic identification, and antimicrobial sensitivity of Aeromonas species isolated from farmed rainbow trout (Onchorynchus mykiss) in Mexico. Acta Trop.; 2014; 130, pp. 76-79. [DOI: https://dx.doi.org/10.1016/j.actatropica.2013.10.021]

12. Zhu, M.; Wang, X.R.; Li, J.; Li, G.Y.; Liu, Z.P.; Mo, Z.L. Identification and virulence properties of Aeromonas veronii bv. sobria isolates causing an ulcerative syndrome of loach Misgurnus anguillicaudatus. J. Fish Dis.; 2016; 39, pp. 777-781. [DOI: https://dx.doi.org/10.1111/jfd.12413] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26525163]

13. Pablos, M.; Huys, G.; Cnockaert, M.; Rodríguez-Calleja, J.M.; Otero, A.; Santos, J.A.; García-López, M.L. Identification and epidemiological relationships of Aeromonas isolates from patients with diarrhea, drinking water and foods. Int. J. Food Microbiol.; 2011; 147, pp. 203-210. [DOI: https://dx.doi.org/10.1016/j.ijfoodmicro.2011.04.006]

14. Chen, P.L.; Tsai, P.J.; Chen, C.S.; Lu, Y.C.; Chen, H.M.; Lee, N.Y.; Lee, C.C.; Li, C.W.; Li, M.C.; Wu, C.J. et al. Aeromonas stool isolates from individuals with or without diarrhea in southern Taiwan: Predominance of Aeromonas Veronii. J. Microbiol. Immunol. Infect.; 2015; 48, pp. 618-624. [DOI: https://dx.doi.org/10.1016/j.jmii.2014.08.007]

15. Parker, J.L.; Shaw, J.G. Aeromonas spp. clinical microbiology and disease. J. Infect.; 2011; 62, pp. 109-118. [DOI: https://dx.doi.org/10.1016/j.jinf.2010.12.003] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21163298]

16. Zhu, X.; Qin, L.; Zhu, Y.; Qian, Q.; Gao, X.; Jiang, Q.; Wang, J.; Liu, G.; Zhang, X. Characteristics and Complete Genome Analysis of a Pathogenic Aeromonas Veronii SJ4 from Diseased Siniperca Chuatsi. Mar. Biotechnol.; 2023; 25, pp. 966-982. [DOI: https://dx.doi.org/10.1007/s10126-023-10253-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37947961]

17. Liu, G.; Li, J.; Jiang, Z.; Zhu, X.; Gao, X.; Jiang, Q.; Wang, J.; Wei, W.; Zhang, X. Pathogenicity of Aeromonas veronii causing mass mortalities of Odontobutis potamophila and its induced host immune response. Fish Shellfish Immunol.; 2022; 125, pp. 180-189. [DOI: https://dx.doi.org/10.1016/j.fsi.2022.05.009] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35561950]

18. Li, Z.; Wang, X.; Chen, C.; Gao, J.; Lv, A. Transcriptome profiles in the spleen of African catfish (Clarias gariepinus) challenged with Aeromonas Veronii. Fish Shellfish Immunol.; 2019; 86, pp. 858-867. [DOI: https://dx.doi.org/10.1016/j.fsi.2018.12.029]

19. Wang, X.; Liu, L.; Zhang, R.; Li, H.; Zhu, H. Involvement of galectin-9 from koi carp (Cyprinus carpio) in the immune response against Aeromonas veronii infection. Fish Shellfish Immunol.; 2022; 129, pp. 64-73. [DOI: https://dx.doi.org/10.1016/j.fsi.2022.08.006]

20. Tang, H.; Zhang, J.; Zhu, L.; Jiang, X.; Pei, C.; Li, L.; Kong, X. Characteristics of CD4-1 gene and its immune responses against Aeromonas veronii infection by activating NF-κB signaling in Qihe crucian carp Carassius Auratus. Fish Shellfish Immunol.; 2024; 145, 109318. [DOI: https://dx.doi.org/10.1016/j.fsi.2023.109318]

21. Tang, H.; Zhu, L.; Zhao, X.; Jiang, X.; Zhang, J.; Pei, C.; Li, L.; Kong, X. Characterization of CD3γ/δ gene and its immune response in Qihe crucian carp Carassius auratus after challenged by Aeromonas veronii and Poly(I:C). Fish Shellfish Immunol.; 2023; 133, 108550. [DOI: https://dx.doi.org/10.1016/j.fsi.2023.108550]

22. Han, Z.; Sun, J.; Jiang, B.; Hu, X.; Lv, A.; Chen, L.; Guo, Y. Concurrent infections of Aeromonas veronii and Vibrio cholerae in koi carp (Cyprinus carpio Var. Koi). Aquaculture; 2021; 535, 736395. [DOI: https://dx.doi.org/10.1016/j.aquaculture.2021.736395]

23. Chen, F.; Sun, J.; Han, Z.; Yang, X.; Xian, J.A.; Lv, A.; Hu, X.; Shi, H. Isolation, Identification and Characteristics of Aeromonas veronii From Diseased Crucian Carp (Carassius auratus gibelio). Front. Microbiol.; 2019; 10, 2742. [DOI: https://dx.doi.org/10.3389/fmicb.2019.02742] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32038507]

24. Hoai, T.D.; Trang, T.T.; Van Tuyen, N.; Giang, N.T.H.; Van Van, K. Aeromonas veronii caused disease and mortality in channel catfish in Vietnam. Aquaculture; 2019; 513, 734425. [DOI: https://dx.doi.org/10.1016/j.aquaculture.2019.734425]

25. Raj, N.S.; Swaminathan, T.R.; Dharmaratnam, A.; Raja, S.A.; Ramraj, D.; Lal, K.K. Aeromonas veronii caused bilateral exophthalmia and mass mortality in cultured Nile tilapia, Oreochromis niloticus (L.) in India. Aquaculture; 2019; 512, 734278. [DOI: https://dx.doi.org/10.1016/j.aquaculture.2019.734278]

26. Zwollo, P.; Cole, S.; Bromage, E.; Kaattari, S. B cell heterogeneity in the teleost kidney: Evidence for a maturation gradient from anterior to posterior kidney. J. Immunol.; 2005; 174, pp. 6608-6616. [DOI: https://dx.doi.org/10.4049/jimmunol.174.11.6608]

27. Bromage, E.S.; Kaattari, I.M.; Zwollo, P.; Kaattari, S.L. Plasmablast and plasma cell production and distribution in trout immune tissues. J. Immunol.; 2004; 173, pp. 7317-7323. [DOI: https://dx.doi.org/10.4049/jimmunol.173.12.7317]

28. Qin, G.; Ai, X.; Xu, J.; Yang, Y. Dual RNA-seq of spleens extracted from channel catfish infected with Aeromonas veronii reveals novel insights into host-pathogen interactions. Ecotoxicol. Environ. Saf.; 2023; 252, 114609. [DOI: https://dx.doi.org/10.1016/j.ecoenv.2023.114609]

29. Han, C.; Li, Q.; Chen, Q.; Zhou, G.; Huang, J.; Zhang, Y. Transcriptome analysis of the spleen provides insight into the immunoregulation of Mastacembelus armatus under Aeromonas veronii infection. Fish Shellfish Immunol.; 2019; 88, pp. 272-283. [DOI: https://dx.doi.org/10.1016/j.fsi.2019.02.020]

30. Choi, K.C.; Lee, Y.S.; Lim, S.; Choi, H.K.; Lee, C.H.; Lee, E.K.; Hong, S.; Kim, I.H.; Kim, S.J.; Park, S.H. Smad6 negatively regulates interleukin 1-receptor-Toll-like receptor signaling through direct interaction with the adaptor Pellino-1. Nat. Immunol.; 2006; 7, pp. 1057-1065. [DOI: https://dx.doi.org/10.1038/ni1383]

31. Burge, S.; Kelly, E.; Lonsdale, D.; Mutowo-Muellenet, P.; McAnulla, C.; Mitchell, A.; Sangrador-Vegas, A.; Yong, S.Y.; Mulder, N.; Hunter, S. Manual GO annotation of predictive protein signatures: The InterPro approach to GO curation. Database; 2012; 2012, bar068. [DOI: https://dx.doi.org/10.1093/database/bar068]

32. Yulis, M.; Quiros, M.; Hilgarth, R.; Parkos, C.A.; Nusrat, A. Intracellular Desmoglein-2 cleavage sensitizes epithelial cells to apoptosis in response to pro-inflammatory cytokines. Cell Death Dis.; 2018; 9, 389. [DOI: https://dx.doi.org/10.1038/s41419-018-0380-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29523777]

33. Fu, M.; Blackshear, P.J. RNA-binding proteins in immune regulation: A focus on CCCH zinc finger proteins. Nat. Rev. Immunol.; 2017; 17, pp. 130-143. [DOI: https://dx.doi.org/10.1038/nri.2016.129] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27990022]

34. Gao, N.; Hibi, Y.; Cueno, M.; Asamitsu, K.; Okamoto, T. A-kinase-interacting protein 1 (AKIP1) acts as a molecular determinant of PKA in NF-κB signaling. J. Biol. Chem.; 2010; 285, pp. 28097-28104. [DOI: https://dx.doi.org/10.1074/jbc.M110.116566] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20562110]

35. Veillette, A.; Dong, Z.; Latour, S. Consequence of the SLAM-SAP signaling pathway in innate-like and conventional lymphocytes. Immunity; 2007; 27, pp. 698-710. [DOI: https://dx.doi.org/10.1016/j.immuni.2007.11.005]

36. Detre, C.; Keszei, M.; Romero, X.; Tsokos, G.C.; Terhorst, C. SLAM family receptors and the SLAM-associated protein (SAP) modulate T cell functions. Semin. Immunopathol.; 2010; 32, pp. 157-171. [DOI: https://dx.doi.org/10.1007/s00281-009-0193-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20146065]

37. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics; 2014; 30, pp. 2114-2120. [DOI: https://dx.doi.org/10.1093/bioinformatics/btu170]

38. Li, J.T.; Wang, Q.; Huang Yang, M.D.; Li, Q.S.; Cui, M.S.; Dong, Z.J.; Wang, H.W.; Yu, J.H.; Zhao, Y.J.; Yang, C.R. et al. Parallel subgenome structure and divergent expression evolution of allo-tetraploid common carp and goldfish. Nat. Genet.; 2021; 53, pp. 1493-1503. [DOI: https://dx.doi.org/10.1038/s41588-021-00933-9]

39. Mortazavi, A.; Williams, B.A.; McCue, K.; Schaeffer, L.; Wold, B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods; 2008; 5, pp. 621-628. [DOI: https://dx.doi.org/10.1038/nmeth.1226]

40. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol.; 2015; 33, pp. 290-295. [DOI: https://dx.doi.org/10.1038/nbt.3122]

41. Mistry, J.; Chuguransky, S.; Williams, L.; Qureshi, M.; Salazar, G.A.; Sonnhammer, E.L.L.; Tosatto, S.C.E.; Paladin, L.; Raj, S.; Richardson, L.J. et al. Pfam: The protein families database in 2021. Nucleic Acids Res.; 2021; 49, pp. D412-D419. [DOI: https://dx.doi.org/10.1093/nar/gkaa913]

42. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res.; 2000; 28, pp. 27-30. [DOI: https://dx.doi.org/10.1093/nar/28.1.27]

43. Trapnell, C.; Williams, B.A.; Pertea, G.; Mortazavi, A.; Kwan, G.; van Baren, M.J.; Salzberg, S.L.; Wold, B.J.; Pachter, L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol.; 2010; 28, pp. 511-515. [DOI: https://dx.doi.org/10.1038/nbt.1621] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20436464]

44. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.; 2014; 15, 550. [DOI: https://dx.doi.org/10.1186/s13059-014-0550-8]

45. Haynes, W. Benjamini–Hochberg Method. Encyclopedia of Systems Biology; Dubitzky, W.; Wolkenhauer, O.; Cho, K.H.; Yokota, H. Springer: New York, NY, USA, 2013; 78.

46. Bu, D.; Luo, H.; Huo, P.; Wang, Z.; Zhang, S.; He, Z.; Wu, Y.; Zhao, L.; Liu, J.; Guo, J. et al. KOBAS-i: Intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis. Nucleic Acids Res.; 2021; 49, pp. W317-W325. [DOI: https://dx.doi.org/10.1093/nar/gkab447]

47. Chen, C.; Wu, Y.; Li, J.; Wang, X.; Zeng, Z.; Xu, J.; Liu, Y.; Feng, J.; Chen, H.; He, Y. et al. TBtools-II: A “one for all, all for one” bioinformatics platform for biological big-data mining. Mol. Plant.; 2023; 16, pp. 1733-1742. [DOI: https://dx.doi.org/10.1016/j.molp.2023.09.010] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37740491]

48. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform.; 2008; 9, 559. [DOI: https://dx.doi.org/10.1186/1471-2105-9-559] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19114008]

49. Zhang, B.; Horvath, S. A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol.; 2005; 4, 17. [DOI: https://dx.doi.org/10.2202/1544-6115.1128]

50. Yip, A.M.; Horvath, S. Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinform.; 2007; 8, 22. [DOI: https://dx.doi.org/10.1186/1471-2105-8-22]

51. Li, A.; Horvath, S. Network neighborhood analysis with the multi-node topological overlap measure. Bioinformatics; 2007; 23, pp. 222-231. [DOI: https://dx.doi.org/10.1093/bioinformatics/btl581]

52. Horvath, S.; Dong, J. Geometric interpretation of gene coexpression network analysis. PLoS Comput. Biol.; 2008; 4, e1000117. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1000117]

53. Xiao, Z.; Shen, J.; Feng, H.; Liu, H.; Wang, Y.; Huang, R.; Guo, Q. Characterization of two thymosins as immune-related genes in common carp (Cyprinus carpio L.). Dev. Comp. Immunol.; 2015; 50, pp. 29-37. [DOI: https://dx.doi.org/10.1016/j.dci.2015.01.003] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25596145]

54. Livak, K.J.; Schmittgen, T.D. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2^−ΔΔCT Method. Methods; 2001; 25, pp. 402-408. [DOI: https://dx.doi.org/10.1006/meth.2001.1262] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11846609]