Dual RNA-seq provides novel insight into the roles of dksA from Pseudomonas plecoglossicida in pathogen-host interactions with large yellow croakers ( Larimichthys crocea).

Pseudomonas plecoglossicida is a rod-shaped, gram-negative bacterium with flagella. It causes visceral white spot disease and high mortality in Larimichthys crocea during culture, resulting in serious economic loss. Analysis of transcriptome and quantitative real-time polymerase chain reaction (PCR) data showed that dksA gene expression was significantly up-regulated after 48 h of infection with Epinephelus coioides (log 2FC=3.12, P<0.001). RNAi of five shRNAs significantly reduced the expression of dksA in P. plecoglossicida, and the optimal silencing efficiency was 96.23%. Compared with wild-type strains, the symptoms of visceral white spot disease in L. crocea infected with RNAi strains were reduced, with time of death delayed by 48 h and mortality reduced by 25%. The dksA silencing led to a substantial down-regulation in cellular component-, flagellum-, and ribosome assembly-related genes in P. plecoglossicida, and the significant up-regulation of fliC may be a way in which virulence is maintained in P. plecoglossicida. The GO and KEGG results showed that RNAi strain infection in L. crocea led to the down-regulation of inflammatory factor genes in immune-related pathways, which were associated with multiple immune response processes. Results also showed that dksA was a virulence gene in P. plecoglossicida. Compared with the wild-type strains, RNAi strain infection induced a weaker immune response in L. crocea.

INTRODUCTION

Infection is an exceedingly complex process involving strong interactions between pathogen and host (Luo et al., 2020). Both pathogen and host must mobilize all available resources to win this life and death battle, and many changes during infection are reflected in their respective transcripts (Zhang et al., 2018a). Therefore, simultaneous detection of transcriptome profiles during infection can provide insight into the pathogenic mechanisms and host immune responses (Tang et al., 2019a). For a long time, due to technical limitations, studies on infection have focused on either the host or the pathogen (Sun et al., 2018). The advancement of dual RNA-seq, which can simultaneously detect both pathogen and host transcriptomes, has provided a powerful and advantageous tool for studying various infection models and pathogen-host interactions (Westermann et al., 2012, 2016, 2017; Valenzuela-Miranda & Gallardo-Escarate, 2018). Dual RNA-seq combined with RNAi has also been used to study the role of virulence genes in host-pathogen interactions (Nuss et al., 2017; Wang et al., 2019a). In addition, dual RNA-seq and dual iTRAQ have been applied to explore gene functions at the multi-omics level (Luo et al., 2019a).

Large yellow croakers (Larimichthys crocea), an economically important marine fish in China, are widely cultured in the Fujian and Zhejiang provinces (Yang et al., 2016). The most important factor threatening L. crocea culture is the frequent occurrence of infection epidemics (Tang et al., 2019c). Visceral white spot disease, one of the most destructive diseases in L. crocea, is caused by Pseudomonas plecoglossicida, a short rod-shaped gram-negative bacterium (Huang et al., 2018; Zhang et al., 2014). Pseudomonas plecoglossicida was first isolated from ayu fish (Plecoglossus altivelis) suffering from bacterial hemorrhagic ascites (Nishimori et al., 2000). In view of its potential to cause great harm to aquaculture, the pathogenic mechanism of P. plecoglossicida has attracted considerable attention (Huang et al., 2019; Tao et al., 2016) and several virulence genes have been identified (Tao et al., 2020; Zhang et al., 2017). Earlier studies reported that P. plecoglossicida commonly circulates in orange-spotted groupers (Epinephelus coioides) and can easily adapt to and proliferate in the splenic environment, resulting in significantly higher pathogen loads in spleens (>200 times) than in other tissues (Luo et al., 2020). In addition, P. plecoglossicida infection can cause irreversible disruption of gut microbiota, resulting in increasing mortality (Li et al., 2020). Furthermore, CspA1 is known to contribute to P. plecoglossicida virulence in a temperature-specific manner via regulation of sigX expression (Huang et al., 2020).

In previous studies from our laboratory, we found that the dksA gene from P. plecoglossicida was highly expressed during host infection (data deposited in GenBank SRA database under accession numbers SRP114910 and SRP115064). The dksA gene encodes an RNA polymerase-binding transcription factor (Kamarthapu et al., 2016). The gene has a variety of important functions, such as regulating rRNA promoter activity (Paul et al., 2004), HFQ gene expression and virulence of Shigella flexneri (Sharma & Payne, 2006), growth rate of Escherichia coli (Mallik et al., 2006), and assembly of flagella (Dalebroux et al., 2010). Thus, based on our previous laboratory research and reported literature, dksA may play a role in the pathogenicity of P. plecoglossicida. To date, however, there are no reports on the function of dksA in pathogen-host interactions.

To explore the roles of dksA in host-pathogen interactions between P. plecoglossicida and L. crocea, a dksA-silenced strain of P. plecoglossicida was constructed by RNAi technology and differences in virulence between wild-type and RNAi strains were analysed. The spleens of L. crocea infected by the wild-type or dksA-RNAi strains of P. plecoglossicida were subjected to dual RNA-seq. The present study provides novel insight into host-pathogen interactions between P. plecoglossicida and L. crocea.

MATERIALS AND METHODS

Bacterial strains and culture conditions

The highly pathogenic wild-type strain of P. plecoglossicida (NZBD9) was isolated from the spleen of a diseased large yellow croaker (Huang et al., 2018). Pseudomonas plecoglossicida was shake-cultured (220 r/min) in Luria Bertani (LB) broth at 18 or 28 °C. In addition, Escherichia coli DH5α was obtained from the Beijing Tiangen Company (China) and cultured in LB broth at 37 °C and 220 r/min.

Construction of P. plecoglossicida RNAi strain

The RNAi strains were constructed according to Sun et al. (2018). Five short hairpin RNA sequences targeting the dksA gene were designed using the RNAi website (http://rnaidesigner.thermofisher.com/rnaiexpress/setOption.do?designOption=shrna&pid=708587103220684543), and then synthesized by Shanghai Generay Biotech Co., Ltd. (China) (Supplementary Table S1). Each oligonucleotide was annealed and ligated to the pCM130/tac vector linearized with the restriction enzymes NsiI and BsrGI (New England Biolabs, USA) using T4 DNA ligase (New England Biolabs, USA) (Guo et al., 2018). The preparation of E. coli DH5α competent cells was performed by the CaCl2 method (Mandel & Higa, 1970). The recombinant pCM130/tac vector was transformed into competent E. coli DH5α by heat shock, and then extracted for electroporating into P. plecoglossicida competent cells (Tang et al., 2019b). Finally, the expression of dksA in five dksA-RNAi strains of P. plecoglossicida was verified by quantitative real-time polymerase chain reaction (qRT-PCR).

Artificial infection and sampling

All fish experiments were carried out strictly following the “Guide for the Care and Use of Laboratory Animals” established by the National Institutes of Health. All animal protocols were approved by the Jimei University Animal Ethics Committee (Acceptance No. JMULAC201159).

In total, 200 healthy L. crocea (average weight ~50 g) were obtained from Ningde (Fujian, China) and acclimatized to laboratory conditions for one week. To detect the pathogenicity of P. plecoglossicida in L. crocea, 60 fish were randomly divided into three groups. Each acclimatized fish was intrapleurally injected with 104 colony forming units per gram fish (cfu/g) of the wild-type or RNAi strain of P. plecoglossicida, respectively. As the negative control, other L. crocea individuals were intrapleurally injected with 0.1 mL of phosphate-buffered saline (PBS). The water temperature throughout the experiment was maintained at 18±1 °C. The status of injected fish was recorded twice a day. For spleen sampling, six spleens from wild-type or RNAi strain-infected L. crocea were randomly sampled at 48 h post infection (hpi) for dual RNA-seq, with two spleens mixed as one sample. In addition, six spleens from L. crocea infected by the wild-type or RNAi strain were randomly sampled at 6, 12, 24, 48, 72, and 96 hpi for pathogen load and dksA expression assays.

DNA and RNA isolation

DNA purification of spleen samples was performed according to the instructions provided with the EasyPure Marine Animal Genomic DNA Kit (TransGen Biotech, China). The extracted genomic DNA was stored at –20 °C until use.

Total RNA was extracted using an Eastep® Super Total RNA Extraction Kit (Shanghai Promega Biological Products, Ltd., China). The quality of total RNA was checked by agarose gel electrophoresis. cDNA was synthesized by TransScript All-in-One First-Strand cDNA Synthesis SuperMix for qPCR (One-Step gDNA Removal) (TransGen Biotech, China) (Liu et al., 2017). The synthesized cDNA was used as a new sample template for qRT-PCR and then stored at –20 °C until use.

qRT-PCR

qRT-PCR was performed using a QuantStudio 6 Flex Real-Time PCR System (Life Technologies, USA). All primer sequences were designed using Primer Premier 5.0 (Supplementary Table S2). Bacterial gene expression was normalized using 16S rDNA, L. crocea gene expression was normalized using β-actin, and the relative level of gene expression was calculated using the 2-ΔΔCt method. To assess the pathogen load of P. plecoglossicida in the infected spleens, the gyrB gene copy number was used to assess the number of P. plecoglossicida.

Dual RNA-seq and transcriptome data analysis

Library preparation and sequencing

Sequencing experiments were carried out using an Illumina TruseqTM RNA Sample Prep Kit (Illumina, USA). Total RNA was extracted from tissue samples using TRIzol® reagent, with concentration and purity measured using a Nanodrop 2000, RNA integrity detected by agarose gel electrophoresis, and RIN values determined by an Agilent 2100. A single database requires a total RNA of 1 μg, concentration of ≥50 ng/μL, and OD260/280 of between 1.8 and 2.2. Eukaryotic and prokaryotic mRNAs were simultaneously obtained by removing rRNA, and fragmentation was carried out by adding a fragmentation buffer. After reverse-synthesizing the cDNA, an end repair mix was added to make it blunt-ended, and then poly(A) was added to the 3' end for ligation of the Y-shaped linker. Sequencing was performed on an Illumina HiSeq4000 sequencing platform from Majorbio Biotech Co., Ltd. (China). The RNA sequence data were deposited in the GenBank SRA database under accession Nos. PRJNA607373 and SRP176599. The basis of transcriptome data analysis is high-quality sequencing. Analysis showed that the distribution of the A/T/G/C base content was uniform (Supplementary Figure S2), the base mass distribution of the sequence data met the requirements of subsequent analysis (Supplementary Table S3), and the repeatability between the three repeated samples was good (Supplementary Figure S1).

Processing and mapping of reads

Trimming and quality control of raw Illumina reads were performed using Sickle (https://github.com/najoshi/sickle) and SeqPrep (https:// github.com/jstjohn/SeqPrep) with default settings. For RNA-seq, clean data were mapped to the genome of P. plecoglossicida strain NB2011 (NCBI RefSeq accession number: NZ_ASJX00000000.1) using Bowtie2 (Langmead & Salzberg, 2012). Clean data were mapped to the genome of L. crocea (NCBI RefSeq accession number: GCF_000972845.2) using Hisat2 (Kim et al., 2015).

Differentially expressed mRNAs (DEMs) and enrichment analysis

Differential expression was determined using edgeR, which performs differential expression calculations based on mRNA read count data and a negative binomial distribution model (Anders & Huber, 2010; Robinson et al., 2010). The screening criteria for significant DEMs were: FDR<0.05 and |log2FC|>=1. The DEMs were then subjected to enrichment analysis by hypergeometric distribution testing using Goatools (https://github.com/tanghaibao/goatools) and KOBAS (http://kobas.cbi.pku.edu. cn/home.do) (Xie et al., 2011).

Statistical analyses

All data are expressed as means±standard deviation (SD) from at least three sets of independent experiments. Data analysis was performed using SPSS 18.0 (SPSS Inc., USA), and one-way analysis of variance with Dunnett’s test was used. P-values of <0.05 were considered statistically significant.

RESULTS

Construction of dksA-RNAi strain

Based on qRT-PCR, the expression level of dksA in the P. plecoglossicida-infected spleens at 48 hpi was six times higher than that in the in vitro culture, consistent with the transcriptome analysis results (Figure 1A). The expression level of dksA was significantly down-regulated in the five mutant strains. The mRNA levels in the shRNA-31, shRNA-49, shRNA-81, shRNA-87, and shRNA-249 mutant strains were only 32.19%, 30.80%, 61.18%, 3.77%, and 9.80% that of the wild-type strain, respectively (Figure 1B). The shRNA-87 mutant strain (hereinafter referred to as the dksA-RNAi or RNAi strain) exhibited the lowest dksA RNA level and was thus selected for subsequent analysis. The growth rates of the dksA-RNAi strain and wild-type strain of P. plecoglossicida were determined, although no significant differences between the two strains were observed (Figure 1C).

Figure 1

Construction and characterization of dksA-RNAi strains of P. plecoglossicida

A: Relative expression of gene dksA (in vivo/in vitro). B: dksA expression level in five mutant strains. C: Growth curve of two strains of P. plecoglossicida.

Effect of dksA on P. plecoglossicida pathogenicity

Compared with the wild-type strain, the dksA-RNAi strain of P. plecoglossicida exhibited a significant decrease in virulence, as observed by the 25% increase in the survival rate of infected L. crocea and 48 h delay in first death. No deaths were recorded in the negative control group of L. crocea injected with PBS (Figure 2A). Furthermore, there were significant differences in spleen appearance between the two groups of L. crocea infected by the dksA-RNAi or wild-type strains of P. plecoglossicida. The spleens of L. crocea infected by the wild-type strain showed a large number of typical white nodules on the surface at 60 hpi, whereas the spleens of L. crocea infected by the dksA-RNAi strain displayed only a small number of white spots on the surface (Figure 2B).

Figure 2

Infection of L. crocea by two strains of P. plecoglossicida

A: Survival rate of L. crocea infected by P. plecoglossicida. B: Spleen appearance of infected L. crocea. C: Relative pathogen load of RNAi strain in spleen of L. crocea compared to that of wild-type strain. D: Expression levels of dksA gene in two strains of P. plecoglossicida in vitro (0 hpi) and in vivo (24, 48, 72 and 96 hpi).

Throughout the infection process, the pathogen load in the spleens of L. crocea infected with the dksA-RNAi strain was always lower than that in L. crocea infected with the wild-type strain of P. plecoglossicida, and showed a tendency to increase gradually with the increase in infection time (Figure 2C). In vivo expression levels of dksA in the dksA-RNAi and wild-type strains of P. plecoglossicida were high throughout the infection process, and the expression levels in the dksA-RNAi strain were always lower than that in the wild-type strain (Figure 2D).

Analysis and verification of transcriptome data of P. plecoglossicida

We used edgeR software to calculate gene expression levels. The screening criteria for significant differentially expressed genes (DEGs) were FDR<0.05 and |log2FC|>=1. From the constructed volcano map, a total of 4 988P. plecoglossicida mRNAs were obtained from the transcriptome of the L. crocea spleens infected with the dksA-RNAi strain. Compared with the wild-type strain, we identified 145 differentially expressed P. plecoglossicida mRNAs in the dksA-RNAi strain-infected spleens, 24 of which were up-regulated and 121 of which were down-regulated (Figure 3A). The heat map (Figure 3B) shows pathogenic genes whose differential expression in the three samples exhibited >2-fold change by transcriptome sequencing. The most variable up-regulated gene was L321_RS17380 (log2FC=5.07), and the gene with the largest down-regulated fold-change was flgC (log2FC=–15.02) (Figure 3B). Five up-regulated and down-regulated DEGs of the pathogen were randomly selected for qRT-PCR detection. The qRT-PCR results were consistent with the transcriptome sequencing results (Figure 3C).

Figure 3

DEG enrichment analysis of pathogen transcriptome data

A: Volcano plot obtained from edgeR analysis of P. plecoglossicida transcriptome. B: Heat map of DEGs from transcripts of pathogen in host spleen (FDR<0.05, |log2FC|≥1), C: Five up-regulated and down-regulated DEGs were randomly selected for qRT-PCR verification.

Enrichment analysis of DEGs of P. plecoglossicida

Gene Ontology (GO) enrichment analysis of P. plecoglossicida was performed based on Goatools using Fisher’s exact test. Results showed that there were 292 enriched terms, including 64 significantly enriched terms related to cellular component (12), biological process (45), and molecular function (7). The top 10 terms were organelle part, cell part, cellular nitrogen compound biosynthetic process, bacterial-type flagellum-dependent cell motility, structural molecule activity, cilium or flagellum-dependent cell motility, archaeal or bacterial-type flagellum-dependent cell motility, cell motility, movement of cell or subcellular component, and cellular process (Figure 4A). KOBAS was applied for KEGG pathway enrichment analysis, with 69 terms found to be enriched, including three significantly enriched terms (i.e., flagellar assembly, photosynthesis, and ribosome) (Figure 4B). In total, DEGS were enriched in 30 terms (P<0.01), as displayed in detail inFigure 4C. Compared with the wild-type strain, most DEGs were down-regulated, except for L321-RS02785, L321-RS02790, L321-RS14705, L321-RS17525, L321-RS12450, and L321-RS11375, which were up-regulated and primarily enriched in single-organism cellular process (GO:0044763) and cellular process (GO:0009987) (Figure 4C).

Figure 4

GO and KEGG enrichment analysis of pathogen transcriptome

A: GO enrichment analysis of pathogen DEGs. B: KEGG pathway enrichment analysis of pathogen DEGs. C: Heat map of DEGs for top 30 significantly enriched GO terms.

Analysis and verification of transcriptome data of infected L. crocea

The software used for gene expression calculation and the screening criteria for significant DEGs were the same as those used for the pathogen. From the constructed volcanic map, 27 520 L. crocea mRNAs were obtained from the spleen transcriptome of L. crocea infected with the dksA-RNAi strain. Compared with the wild-type-infected L. crocea, we identified 970 differentially expressed P. plecoglossicida mRNAs in the dksA-RNAi strain-infected spleens, 366 of which were up-regulated and 604 of which were down-regulated (Figure 5A). The heat map (Figure 5B) shows the expression levels of the top 50 up-regulated and down-regulated genes in the three samples of each group. The greatest fold-change in the up-regulated genes was for LOC104939330 (log2FC=8.51), and greatest fold-change in the down-regulated genes was for CPA2 (log2FC=–5.4). Five up-regulated and down-regulated DEGs of the host were randomly selected for qRT-PCR detection. The qRT-PCR results of these genes were consistent with the transcriptome sequencing results (Figure 5C).

Figure 5

DEG enrichment analysis of L. crocea transcriptome data

A: Volcano plot obtained from edgeR analysis of infected L. crocea RNA pools. B: Heat map analysis of spleen transcript DEGs in infected L. crocea (FDR<0.05, |log2FC|≥1), C: Five up-regulated and down-regulated DEGs of host were randomly selected for qRT-PCR verification.

Enrichment analysis of DEGs in L. crocea

GO enrichment analysis of DEGs in L. crocea identified 292 enriched terms, including 21 significantly enriched terms related to cellular component (1), biological process (9), and molecular function (11) (Figure 6A). The top 10 significantly enriched terms for the DEGs are displayed in detail in Figure 6B. Based on z-score, more down-regulated genes were enriched in endopeptidase activity (GO: 0004175), with peptidase activity acting on L-amino acid peptides (GO: 0070011) showing the highest significance (Figure 6B). KEGG analysis identified 265 enriched KEGG pathways, including four significantly enriched immune-related pathways (i.e., Toll-like receptor signaling pathway, tumor necrosis factor (TNF) signaling pathway, hematopoietic cell lineage, and cytokine-cytokine receptor interaction) (Figure 6C). According to the heat map, compared with the wild-type strain-infected spleens, 48 DEGs (39 down-regulated, nine up-regulated) were enriched in the cytokine-cytokine receptor interaction (ko04060) pathway. We identified 14 DEGs (11 down-regulated, three up-regulated) enriched in the hematopoietic cell lineage (ko04640) pathway. In addition, 17 DEGs (16 down-regulated genes, one up-regulated) were enriched in the Toll-like receptor signaling pathway (ko04620) (Figure 8D).

Figure 6

GO and KEGG enrichment analysis of host transcriptome

A: GO enrichment analysis of pathogen DEGs. B: Top 10 GO term-gene annotation enrichment analysis results. C: KEGG pathway enrichment analysis of host DEGs. D: DEGs of cytokine-cytokine receptor interaction (Ko04060), hematopoietic cell lineage (Ko04640). and Toll-like receptor signaling pathway (Ko04620).

According to the mapped TNF signaling pathway, many genes in the spleen infected with the silenced strain showed significant changes in pathways compared to the spleen infected with the wild-type strain, with 23 genes found to be significantly down-regulated. Intracellular signaling (negative), transcription factors, leukocyte recruitment, cell adhesion, and other related gene expression levels changed, with major changes in inflammatory cytokine genes IL1b (–5.22), IL6 (–4.02) and synthesis of inflammatory mediator gene ptgs2 (–4.52) (Figure 7).

Figure 7

Schematic overview of response of TNF signaling pathway of L. crocea to P. plecoglossicida infection

Color indicates relative expression level of gene (RNAi strain/wild-type strain), red is up-regulated, blue is down-regulated. Darker color indicates greater change in gene expression.

DISCUSSION

The ability of pathogens to infect hosts is mainly regulated by virulence genes (Crofts et al., 2018; Yao et al., 2019). In recent years, dozens of virulence genes of aquatic pathogens have been identified (Rong et al., 2017; Zhang et al., 2018b), including several virulence genes of P. plecoglossicida (Wang et al., 2019b). Virulence genes can affect host immune responses (Sun et al., 2018) and are involved in host-pathogen interactions (Sun et al., 2019a; Tang et al., 2020). To date, however, no studies have reported on the effects of the dksA gene on host immune response.

RNA interference technology can specifically reduce the expression of genes and has been widely used to explore gene functions (Zhang et al., 2019a). RNA interference technology results in different silencing efficiency of different aquatic pathogens (Guanzon & Maningas, 2018; Saleh et al., 2016), and different shRNAs have different silencing efficiency for the same gene (Sun et al., 2019b). To achieve good silencing, it is necessary to design several different shRNAs for a gene. In the present study, five shRNAs were designed to silence the dksA gene by RNAi. Among the five RNAi strains, dksA-shRNA-87 had the highest silencing efficiency (96.23%), which is more efficient than the silencing of most genes of aquatic pathogens (Ye et al., 2018; Zuo et al., 2019). The stability of gene silencing is crucial for the study of gene function. In the present study, the expression of dksA in the host was higher than that in vitro, indicating that the gene may be related to the pathogenicity of P. plecoglossicida. Moreover, the dksA gene in the dksA-RNAi strain of P. plecoglossicida was persistently silenced during the infection process and its relative expression was always lower than that of the wild-type strain. These results indicate that the RNAi technique was reliable in this study, thereby laying the foundation for subsequent research.

In this study, L. crocea was artificially infected with P. plecoglossicida. After infection with the dksA-silenced strain, death time was delayed by 48 h and the mortality rate was reduced by 25% compared with L. crocea infected with wild-type P. plecoglossicida. After 60 h of infection, the spleens of L. crocea showed typical symptoms of visceral white spot disease, with fewer symptoms observed in the RNAi strain. RNAi of the dksA gene resulted in a decrease in the pathogenicity of P. plecoglossicida to L. crocea, suggesting that dksA may contribute to the virulence of P. plecoglossicida. Several genes have been verified to be associated with P. plecoglossicida virulence, the silencing of which results in a reduction in mortality of experimental fish (Tang et al., 2019a). The silencing of some genes also results in a dramatic decline in the virulence of pathogens, and therefore shows potential in the development of attenuated vaccines (Luo et al., 2020).

Transcriptome sequencing is a powerful tool for studying gene function (Yang et al., 2018), and has been applied in host-pathogen interaction studies (Zhang et al., 2019b). In this study, the transcriptome of the silenced strain changed significantly during the infection process, and a total of 4 988 P. plecoglossicida mRNAs were detected in the samples. Compared with the wild-type strain, we identified 145 differentially expressed mRNAs in the dksA-RNAi strain-infected spleens, 24 of which were up-regulated and 121 of which were down-regulated. Significantly altered genes were analyzed for GO and KEGG enrichment. As a result, cellular component-, flagellum-, and ribosome assembly-related genes of P. plecoglossicida were down-regulated after dksA gene silencing, and the up-regulated gene L321-RS14705 (fliC) (log2FC=3.24) was enriched in seven functional GO terms. Research has shown that DksA and ppGpp in Escherichia coli inhibit the expression of the flagellar cascade during the stationary phase and following starvation, thus affecting flagella and ribosome assembly (Lemke et al., 2009). The loss of the fliC gene in Edwardsiella tarda can damage bacterial growth, reduce motility, decrease biofilm formation, and decrease secretion of virulence-related proteins involved in the type III secretion system (TTSS) (He et al., 2012). From this perspective, the silencing of the dksA gene in P. plecoglossicida inhibited the expression of a large number of genes related to flagella and ribosome assembly, and the significant up-regulation of the fliC gene may be a way to maintain virulence in P. plecoglossicida.

Infected spleens were chosen for dual RNA-seq because they are an important immune organ (Chen et al., 2019). Dual RNA-seq can synchronously detect transcriptome changes in both host and pathogen (Luo et al., 2019a). During infection of L. crocea, the transcriptome changed significantly. Compared with the wild-type infected group, the most significantly down-regulated gene in the dksA-RNAi strain-infected group was CPA2. CPA2 is a member of the carboxypeptidase gene family. In bacteria, carboxypeptidases play a key role in the immune response to viral infections (Gardell et al., 1988; Godahewa et al., 2014). Based on enrichment analysis, most significantly enriched GO terms were related to peptidase activity, with the greatest impact on endopeptidase activity (GO: 0004175) and peptidase activity acting on L-amino acid peptides (GO: 0070011). Many peptidases are related to the immune response. Endopeptidases and L-amino acid peptides have immune-related functions in organisms, e.g., the role of endopeptidases in the immune response against influenza in mice (Tan et al., 2017) and L-amino acid antibacterial activity in the mucus layer of flounder Platichthys stellatus (Kasai et al., 2010). Some up-regulated genes in the pathogenic bacteria found here have not been reported previously in relation to pathogenicity, and further studies are needed to determine whether they are new virulence genes.

KEGG enrichment analysis of the transcriptome data of L. crocea identified four significantly changed immune-related pathways, i.e., Toll-like receptor signaling pathway, TNF signaling pathway, hematopoietic cell lineage, and cytokine-cytokine receptor interaction. Hematopoietic cell lineage plays an important role in the immune response (Delves, 2020). Through hematopoietic cell lineage, hematopoietic stem cells differentiate into different blood cells, including T cells, natural killer (NK) cells, basophils, macrophages, and B cells, in response to various stimuli (Lu & Chen, 2019; Marshall et al., 2018). The cytokine-cytokine receptor interaction pathway is mainly involved in neutrophil infiltration during host immune response (Mantovani et al., 2019). Cytokines are soluble proteins secreted by donor cells in response to stimuli and transported to target cells through the circulatory system (Sharma et al., 2014). Studies have shown that Epinephelus coioides infected with L321_RS19110 gene-silenced P. plecoglossicida strains can significantly affect cytokine-cytokine receptor interactions (Zhang et al., 2018a). Toll-like receptors play an important role in a host’s ability to recognize pathogens and generate an immune response, which they regulate by promoting inflammatory cytokines (Palti, 2011; Takeda & Akira, 2015). The Toll-like receptor signaling pathway is also involved in the immune response of L. crocea to the secY gene of P. plecoglossicida (Wang et al., 2019a). TNF is a proinflammatory cytokine, which mediates inflammatory responses and regulates immune functions, with abnormal TNF signal transduction also related to inflammatory diseases (Joosten et al., 2016). To ensure body health, an organism will clear senescent and diseased cells by inducing apoptosis through the tumor necrosis factor superfamily (TNFSF) ligand (Collette et al., 2003). Studies have shown that interactions between tumor necrosis factor receptor type 1 (TNFR1) and nuclear factor kappa-B (NF-κB) are essential for maintenance of the TNFR1 pathway activity and activation of inflammatory cytokines that induce leukocyte recruitment (Alcamo et al., 2001). In addition, the TNFR1 pathway can also prepare antibodies for further pathogen clearance (Stokes et al., 2015). In the current study, significantly enriched inflammatory cytokines and cellular regulatory genes were down-regulated in the TNF pathway, similar to results reported in our previous study (Tang et al., 2019c). To summarize, the down-regulation of pro-inflammatory genes involved in the four significantly changed immune-related pathways may be due to the virulence of the RNAi strain being weaker than that of the wild-type strain. Thus, the immune response of the large yellow croakers infected with the RNAi strain was weaker than those infected with the wild-type strain.

CONCLUSIONS

In conclusion, dksA is a virulence gene of P. plecoglossicida. The silencing of dksA resulted in the down-regulation of cell component-, flagellum-, and ribosome assembly-related genes of P. plecoglossicida, thereby reducing the virulence of P. plecoglossicida. Through analysis of transcriptome data, we found that the fliC gene of the RNAi strain was significantly up-regulated in the course of infection, which may be a way in which to maintain the virulence of P. plecoglossicida. In addition, compared with those infected with the wild-type strain, the immune response of L. crocea infected by the RNAi strain was weaker. Most down-regulated GO terms in L. crocea infected with the RNAi strain were related to peptidase activity. KEGG enrichment analysis showed that genes related to inflammatory factors in four immune-related pathways were down-regulated in L. crocea infected with the RNAi strain. Therefore, L. crocea appears to be more resistant to infection by RNAi strains.

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

Q.P.Y., L.X.H., and Y.Q.S. designed the study. Q.P.Y. supervised the analyses. L.Y.W, L.M.Z., and W.Q.Z. performed the fish experiments. Z.X.L. cultivated the bacteria. L.Y.W. extracted DNA and RNA. L.Y.W., Z.X.L., L.X.H., Y.X.Q., and F.W. performed bioinformatics analysis. L.Y.W. wrote the manuscript with input from other authors. Q.P.Y., L.X.H., and Y.Q.S. revised the manuscript. All authors read and approved the final version of the manuscript.