Novel insights into host-pathogen interactions of large yellow croakers ( Larimichthys crocea) and pathogenic bacterium Pseudomonas plecoglossicida using time-resolved dual RNA-seq of infected spleens.

Host-pathogen interactions are highly complex, involving large dynamic changes in gene expression during infection. These interactions are fundamental to understanding anti-infection immunity of hosts, as well as the pathogenesis of pathogens. For bacterial pathogens interacting with animal hosts, time-resolved dual RNA-seq of infected tissue is difficult to perform due to low pathogen load in infected tissue. In this study, an acute infection model of Larimichthys crocea infected by Pseudomonas plecoglossicida was established. The spleens of infected fish exhibited typical symptoms, with a maximum bacterial load at two days post-injection (dpi). Time-resolved dual RNA-seq of infected spleens was successfully applied to study host-pathogen interactions between L. crocea and P. plecoglossicida. The spleens of infected L. crocea were subjected to dual RNA-seq, and transcriptome data were compared with those of noninfected spleens or in vitro cultured bacteria. Results showed that pathogen-host interactions were highly dynamically regulated, with corresponding fluctuations in host and pathogen transcriptomes during infection. The expression levels of many immunogenes involved in cytokine-cytokine receptor, Toll-like receptor signaling, and other immune-related pathways were significantly up-regulated during the infection period. Furthermore, metabolic processes and the use of oxygen in L. crocea were strongly affected by P. plecoglossicida infection. The WGCNA results showed that the metabolic process was strongly related to the entire immune process. For P. plecoglossicida, the expression levels of motility-related genes and flagellum assembly-related genes were significantly up-regulated. The results of this study may help to elucidate the interactions between L. crocea and P. plecoglossicida.

INTRODUCTION

Large yellow croakers (Larimichthys crocea), a marine fish with the largest annual yield in China, are widely cultured in Fujian and Zhejiang provinces (Sun et al., 2019a; Yang et al., 2016). The most significant factor affecting the production of L. crocea is the frequent outbreak of diseases caused by pathogens (Tang et al., 2019a). The “visceral white spot disease” caused by Pseudomonas plecoglossicida is one of the most destructive diseases of L. crocea, causing considerable economic losses (Huang et al., 2018; Zhang et al., 2014). This pathogen not only infects large yellow croakers but is the causative agent of bacterial hemorrhagic ascites disease in Plecoglossus altivelis (Nishimori et al., 2000). In a previous study, P. plecoglossicida infection (104 colony forming units per gram (cfu/g)) resulted in 100% L. crocea death rate at 4.5 days post-injection (dpi) (Tang et al., 2019c), indicating that P. plecoglossicida is more virulent than most other aquatic pathogens, including Vibrio alginolyticus (Huang et al., 2017; Liu et al., 2017) and Aeromonas hydrophila (Zhang et al., 2018b, 2019c). In view of the considerable damage caused by P. plecoglossicida to aquaculture, its pathogenic mechanism has attracted the attention of researchers (Huang et al., 2018; Tao et al., 2016), resulting in the identification of several virulence genes (Luo et al., 2019a; Zhang et al., 2019a).

Infection is a life-and-death struggle between host and pathogen (Zhang et al., 2018a), in which both pathogen and host must mobilize all available resources to win (Luo et al., 2020). Many of the changes that occur in pathogens and hosts during infection will be reflected in their respective transcriptome profiles (Jenner & Young, 2005; Sorek & Cossart, 2010). Traditional RNA-seq has been widely used (Sun et al., 2018; Yang et al., 2018) but cannot detect both host and pathogen transcriptomes simultaneously. However, dual RNA-seq provides analysis of both host and pathogen transcriptome interactions, providing valuable insights into infectious biology (Westermann et al., 2012, 2017). Dual RNA-seq has been used in many bacterial infection studies, including cell-infected and animal-infected models, and has identified many key genes involved in pathogen-host interactions (Aprianto et al., 2016; Minhas et al., 2019; Westermann et al., 2016). For instance, analysis of host and pathogen transcriptional profiles during Staphylococcus aureus infection of two mouse strains demonstrated that the expression of virulence factors is dependent on encountered host resistance (Thänert et al., 2017). Furthermore, using a murine model of pneumococcal infection, Ritchie & Evans (2019) identified key host and bacterial responses to invasion of the pleural space. Dual RNA-seq has also been used to study virulence genes in host-pathogen interactions. (Sun et al., 2019a; Wang et al., 2019). Furthermore, recent study integrated dual RNA-seq and dual iTRAQ (Luo et al., 2019b). Host-pathogen interactions are a complex dynamic process. To date, however, most dual RNA-seq studies of infected tissue have been obtained through static dual RNA-seq (Luo et al., 2020).

To better understand the interactions of L. crocea and P. plecoglossicida, an infection model was established, and the spleens at different stages of infection were subjected to dual RNA-seq to simultaneously monitor the transcriptomes of both L. crocea and P. plecoglossicida. The results of this study should help elucidate host-pathogen interactions between L. crocea and P. plecoglossicida.

MATERIALS AND METHODS

Experimental fish and temporary culture conditions

All healthy and size-matched L. crocea were provided by the State Key Laboratory of Large Yellow Croaker Breeding, Ningde, China. The fish were acclimatized at 18±1 °C for one week under specific pathogen-free laboratory conditions before infection.

Bacterial strain and culture conditions

The pathogenic P. plecoglossicida (NZBD9) strain was isolated from the spleen of naturally infected L. crocea (Huang et al., 2018). The bacteria were routinely overnight cultured in Luria Bertani (LB) broth at 18 or 28 °C with shaking at 220 r/min.

Infection test

All fish experiments were carried out in strict accordance with 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).

For survival assays, acclimatized L. crocea were injected intrapleurally with 104 cfu/g of P. plecoglossicida (NZBD9) (Tang et al., 2019c). The L. crocea fish intrapleurally injected with 0.1 mL of phosphate-buffered saline (PBS) were used as a negative control. The water temperature during the entire experiment was maintained at 18 °C. The mortality of infected L. crocea was recorded twice a day. Dead L. crocea were removed in a timely manner, and their spleens were dissected and photographed.

For dual RNA-seq assays, the spleens of six infected L. crocea were sampled at 1, 2, 3, and 4 dpi. Two spleens were mixed as one sample. PBS-treated L. crocea and in vitro-cultured P. plecoglossicida (18 °C) were used as controls.

RNA isolation, cDNA library construction, and sequencing

Total RNA of spleens and P. plecoglossicida strains was extracted separately using TRIzol reagent (Invitrogen, USA). Turbo DNA-free DNase (Ambion, USA) treatment was performed to remove possible contaminant genomic DNA. The RNA samples were further purified using a Ribo-Zero rRNA Removal Kit (human/mouse/rat or gram-negative bacteria) (Epicentre, USA). The mRNA samples were then fragmented in fragmentation buffer, and a SuperScript double-stranded cDNA synthesis kit (Invitrogen, USA) was used to synthesize cDNA. After modification and purification, library fragments were enriched by polymerase chain reaction (PCR) amplification. An Agilent 2100 Bioanalyzer (Agilent Technologies, USA) was used to detect library quality. The Illumina NovaSeq sequencing platform was used to carry out second-generation sequencing.

Transcriptome data analysis

Sequence data processing and differential gene expression analysis

The measurement of differential gene expression was performed individually for the P. plecoglossicida and L. crocea libraries. Read mapping and transcript abundance quantification (FPKM-normalized expression value) were performed using Bowtie2 (Kim et al., 2015) and RSEM (Li & Dewey, 2011). Differentially expressed mRNAs (DEMs) at different time points were tested using the package edgeR with FDR<0.05 & |log2 fold-change|≥1. ClusterProfiler and Goatools were used for Gene Ontology (GO) enrichment analysis of L. crocea and P. plecoglossicida, respectively, and only categories with P<0.05 were considered to be enriched in the network (Klopfenstein et al., 2018; Yu et al., 2012). The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used for pathway enrichment analysis (Kanehisa & Goto, 2000).

STEM analysis

To analyze changes in gene expression in different infection time samples, Short Time-series Expression Miner (STEM) (v1.3.11) was used to analyze gene expression patterns in L. crocea and P. plecoglossicida (Ernst & Bar-Joseph, 2006). As recommended by the software, the median FPKM in the biological repeat group was selected as the final expression amount. Clustered profiles with P<0.05 were considered significantly expressed, and DEMs were clustered into different expression patterns for subsequent analysis.

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) was performed on L. crocea and P. plecoglossicida data using the WGCNA package in R (Langfelder & Horvath, 2008). WGCNA can help identify functional pathways when studying species interactions (Grote et al., 2017). Ranked by variance from large to small, we chose the top 25% of genes for WGCNA. A network was created based on the results of the pickSoftThreshold function. Hierarchical clustering and dynamic branch cutting were used to identify stable modules of densely interconnected genes. Hub genes were defined by module connectivity, as measured by the absolute values of Pearson’s correlations. Hub genes in a module are considered to have significant functional significance. In this study, key modules were visualized by Cytoscape (http://www.cytoscape.org/).

Quantitative real-time PCR (qRT-PCR)

To confirm the accuracy of data obtained with dual RNA-Seq, some DEMs from bacteria and fish were randomly selected for qRT-PCR. All primers used are listed in Supplementary Table S2. Primers were designed by Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA). qRT-PCR was performed using QuantStudio 6 Flex (Life Technologies, USA). Three replicates were included for each sample, and the β-actin gene of L. crocea served as an internal control to normalize the expression levels; for bacteria, 16S was used for normalization (Wang et al., 2019). The 2−ΔΔCt method was used to calculate the relative level of differential gene expression (Livak & Schmittgen, 2001). All data represent the mean±standard deviation (SD) of three replicates.

Bacterial load determination

DNA purification from spleens was accomplished with an EasyPure Marine Animal Genomic DNA Kit (TransGen Biotech, China) following the manufacturer’s instructions. The abundance of P. plecoglossicida was estimated by qRT-PCR using the copy number of the gyrB gene (Huang et al., 2019). Three replicates were conducted for each experiment.

Statistical analyses

Data are expressed as means±SD. Data analysis was performed with SPSS 17.0 (SPSS Inc., USA). Independent sample t-test and one-way analysis of variance (ANOVA) with Dunnett’s test were used. P-values less than 0.05 were considered significant.

Data access

The RNA sequencing read data were deposited in the GenBank SRA database under accession Nos. SRP176599 and PRJNA597434.

RESULTS

Pathogen load in spleen after infection

As the causative agent of “visceral white spot disease” in L. crocea, P. plecoglossicida exhibited strong pathogenicity to L. crocea. After intrapleural injection of a relatively low dose (104 cfu/g) of P. plecoglossicida, experimental fish died shortly thereafter. The first death of L. crocea was recorded at 2.5 dpi, and the last death was recorded at 4.5 dpi (Tang et al., 2019c). All L. crocea counterparts from the CG that received intrapleural injections of PBS survived the experimental period.

At 3 dpi, the spleens in L. crocea injected with P. plecoglossicida exhibited typical symptoms (spleen surface covered with numerous white spots); however, no white spots were observed on the surface of spleens in the L. crocea CG (Figure 1A).

1

Pathology and pathogen load in spleen of L. crocea infected by P. plecoglossicida NZBD9

A: Symptoms of spleens. Infection group was injected with P. plecoglossicida, and control group (CG) was injected with PBS. B: Pathogen load of P. plecoglossicida in infected L. crocea.

The qRT-PCR results showed that the P. plecoglossicida load in the infected spleens varied with time of infection. At the beginning of infection, the bacterial load was notably low, then increased rapidly, peaked at 2 dpi, and decreased with time (Figure 1B).

Transcriptome of infected spleens reflects different temporal dynamics during infection

An overview of sequencing and mapping results of L. crocea and P. plecoglossicida in infected spleens is provided in Supplementary Table S1. Principal component analysis (PCA) of fish samples revealed a major difference between infected and uninfected samples (Figure 2A). Fish transcript profiles from infected samples could be further assigned into two groups: WT1d and WT2d–Wt3d–WT4d (13.17% of variance; Figure 2A). Similarly, PCA analysis of P. plecoglossicida transcript profiles identified two groups (Figure 2B). Some host and pathogen genes were randomly chosen for qRT-PCR, with results demonstrating the reliability of dual RNA-seq data. Fold-changes obtained by qRT-PCR and dual RNA-seq showed a relatively high correlation for both species: R2=0.7561 for host transcripts and R2=0.749 for pathogenic transcripts (Figure 2C), validating the reliability of the dual RNA-seq data. The main expression of the transcriptomes of the two species is shown in Figure 2E, D.

2

Transcriptional profiles of L. crocea and P. plecoglossicida in infected spleens

A: PCA clustering of full transcriptional profiles of L. crocea. B: PCA clustering of full transcriptional profiles of P. plecoglossicida. C: Validation of RNA-seq data by qRT-PCR. Blue indicates host, orange indicates pathogen, and R-square is marked. D: Hierarchical clustering of major classes of differentially expressed mRNAs (DEMs, P<0.01, log2 fold-change (Log2FC)≥2) in L. crocea transcriptomes. E: Hierarchical clustering of major classes of DEMs (P<0.01, log2 fold-change (Log2FC)≥2) in P. plecoglossicida transcriptomes. Relative expression levels of each transcript (rows) in each sample (column) are shown. FPKM was log2-transformed and then median-centered by transcript. CG: Control group; WT1d: Infected group sampled at 1 dpi; WT2d: Infected group sampled at 2 dpi, and so on. Same below.

Infection-induced changes in gene expression of L. crocea

To analyze the dynamic changes in the host transcriptome during infection, we performed multiple transcriptome comparisons and subsequent analysis. The number of DEMs between any two transcriptomes (CG, WT1d, WT2d, WT3d, WT4d) of L. crocea is shown in Figure 3A. The comparison group with the largest difference in each column was the CG and infected group.

3

Differentially expressed mRNA (DEM) analysis of L. crocea infected by P. plecoglossicida

A: Numbers of DEMs between time points in host. Most different time point across all comparisons (i.e., within a column) is shown in red, most similar time point across all comparisons (i.e., within a column) is shown in blue. B: Heat map of GO enrichment analyses. Shade of color represents size of P-adjusted value, “NA” indicates not significantly enriched in this comparison group. C: KEGG enrichment analysis. Size of circle represents size of enriched gene number, each color represents a pathway.

There was a clear difference in the number of DEMs between the WT1d, WT2d, WT3d and WT4d groups compared with the control. Therefore, we used GO and KEGG enrichment analysis to determine the biological function of the identified DEMs. The GO enrichment analysis results in biological process (BP) and molecular function (MF) are shown in Figure 3B. Compared with the CG, the GO terms enriched in all infected groups were peptide metabolic process, translation, peptide biosynthetic process, protein heterodimerization activity, structural constituent of ribosome, and protein dimerization activity. These terms are mainly related to protein synthesis. An immunologically related term (GO: 0006955: immune response) could be seen after 3 dpi. Notably, many GO terms related to nucleic acid synthesis were enriched only on the second day (Figure 3B). However, on the second day, many fish deaths occurred after infection.

The KEGG enrichment results showed that the ribosome pathway was significantly enriched in the WT1–4d groups (Figure 3C). The infection-related pathway (Herpes simplex virus 1 infection) was enriched, except in WT1d, whereas WT3d and WT4d were enriched in pathways related to host defense, such as the Toll-like receptor signal pathway and cytokine-cytokine receptor interaction.

Infection is a dynamic process. To use time series information effectively, gene expression pattern clustering was applied to mine the expression rule of genes. In the STEM profile, some specific expression patterns appeared due to the influence of biological laws, but some were random. To remove this, we calculated the P-value of each profile through a specific statistical model, with P<0.001 considered to be an authentic and useful research profile.

A total of 50 profiles were aggregated, but only 11 highly significant expression profiles were selected for subsequent analysis (P<0.001), as shown in Supplementary Figure S1A. We used the same color to represent similar expression patterns (Figure 4). Compared with the CG, the profiles of 18-LC, 41-LC, and 43-LC, in which mRNAs were largely up-regulated, are shown in red. The mRNAs in green showed overall down-regulation, but with an up-regulation trend. The mRNAs in purple and yellow modules were generally down-regulated but fluctuated, whereas the mRNAs in the orange profile were generally up-regulated and stable.

4

Trend analysis of DEM enrichment across infectious stages of L. crocea

Ten gene sets clustered by STEM software. Control group (CG) was represented by 0 d, and top left corner shows number of genes in profile.

To determine the biological function of DEMs in each profile, GO classification and KEGG pathway analysis were carried out. From the GO MF enrichment results in Figure 5A, no overlap was observed in GO terms for any profile. Different expression patterns corresponded to different biological processes. The DEMs in profile 0-LC were enriched with many oxygen-related GO terms, whereas profile 38-LC was enriched with many protein-related terms. For the genes in other profiles, a large number were dominant in enzyme activity, binding, and energy.

5

STEM profile enrichment analysis of L. crocea

A: Top five GO molecular function (MF) terms for six profiles. Size of circle represents number of enriched genes, each color indicates significance. B: KEGG pathway enrichment analysis. Different colors are used to represent different profiles, and height of bar graph indicates significance.

KEGG pathway analysis showed that the DEMs in eight profiles were mapped with 19 pathways, as shown in Figure 5B, e.g., profile 0-LC, profile 1-LC, profile 2-LC, and profile 3-LC. There were multiple profiles enriched with ribosomes and spliceosomes. Notably, profile 43-LC was enriched in immune-related pathway cytokine-cytokine receptor interaction.

Host adaptation of P. plecoglossicida to L. crocea

Figure 6A shows the number of differentially expressed genes between any two transcriptomes (CG, WT1d, WT2d, WT3d, WT4d) of P. plecoglossicida. Compared with the CG, WT3d had the greatest number of DEMs in the infected group, reaching 2 115.

6

Differentially expressed mRNA (DEM) analysis of P. plecoglossicida

A: Numbers of DEMs between time points in host. Most different time point across all comparisons (i.e., within a column) is shown in red, most similar time point across all comparisons (i.e., within a column) is shown in blue. B: GO molecular function (MF) enrichment analyses. Different colors represent different groups, and length of bar graph indicates significance. C: KEGG enrichment analyses.

Compared with the CG, the results of GO enrichment are shown in Figure 6B. The CG/WT1d group was dominated by catalytic activity, the CG/WT2d group was related to motor activity, and the CG/WT3d group was related to RNA binding.

For KEGG enrichment, only the CG/WT2d and CG/WT3d groups had results, and they were all related to flagellar assembly and bacterial chemotaxis (Figure 6C).

To analyze the expression mode of DEMs, which played key roles in the infection process, we performed trend analysis of DEMs from the five transcriptomes. A total of 50 profiles were aggregated, with nine highly significant expression profiles then selected for subsequent analysis (P<0.001), as shown inFigure 7 and Supplementary Figure S1. Compared with the CG, the expression levels of genes in modules 18-PP, 43-PP, and 41-PP were generally up-regulated, whereas gene expression levels in modules 3-PP, 1-PP, 0-PP, and 8-PP were down-regulated. For genes in modules 14-PP and 5-PP, gene expression was first down-regulated and then up-regulated.

7

STEM analysis for P. plecoglossicida

A: Six major profiles of gene expression. Top left-hand corner indicates number of DEMs belonging to profile. Lower left-hand corner shows P-value of profile. Colored lines represent DEMs. X-axis represents days after P. plecoglossicida inoculation (dpi). Y-axes represent log2-fold-change in gene expression between treatments. B: GO enrichment analysis of profile 43-PP, 5-PP, 3-PP, and 1-PP. C: Top ten GO terms of profile 18-PP. Size of circle represents quantity of enriched genes, and each color indicates significance. D: KEGG pathway enrichment analysis. Size of circle represents size of enriched gene number, each color represents different groups.

The dynamic process of P. plecoglossicida during infection was elucidated by GO enrichment analysis of expression profiles. Some profiles did not show useful results; therefore, only 18-PP, 43-PP, 41-PP, 5-PP, 3-PP, and 1-PP (Figure 7A) were selected for analysis. The moudle which enriched the most GO terms was profile 18-PP. Compared with the CG, gene expression decreased on the second day, but was up-regulated thereafter. The GO terms enriched in profile 18-PP were mainly related to ribosomes, transcriptional translation, or energy. The top 10 GO terms are shown in Figure 7C. These GO terms were both necessary and fundamental to pathogen activity.

As shown in Figure 7B, only the enriched GO terms in module 43-PP were related to bacterial movement and flagella: i.e., bacterial-type flagellum-dependent cell motility; cilium or flagellum-dependent cell motility; archaeal or bacterial-type flagellum-dependent cell; motility; cell motility; movement of cell or subcellular component; and bacterial-type flagellum part and locomotion. Based on STEM profiles, the general expression of a gene can be observed. Compared with the CG, the genes in profile 43-PP showed an upward trend until the fourth day, after which gene expression began to decrease slightly when most of the hosts died (Figure 7A). Thus, we speculate that flagella were important for P. plecoglossicida throughout the infection period.

The GO enrichment results of profiles 3-PP and 1-PP are shown in Figure 7B. The expression patterns of these two profiles were similar, and showed enrichment in resistance-related GO terms: i.e., response to stress and response to stimulus.

Co-expression of host and pathogen genes

To determine which genes were co-expressed between L. crocea and P. plecoglossicida during infection, we built a co-expression network for the two organisms using WGCNA (Figure 8). In the same module, co-expressed genes at different time points were identified. In total, the CG and four experimental groups contained 15 samples, with variance then calculated. The top 25% of genes were used to construct the gene co-expression network. The choice of soft-thresholding power was an important step in building a co-expression network. We performed network topology analysis of the soft-thresholding power from 1 to 20 and determined the relative balance of scale independence and average connectivity of the network. The middle results are shown in Supplementary Figure S2.

8

Weighted Gene Co-expression Network Analysis (WGCNA) of transcriptional changes in L. crocea and P. plecoglossicida

A: L. crocea time-resolved transcriptome WGCNA cluster map. B: KEGG enrichment analysis for yellow module. Size of circle represents number of genes in pathway, color indicates significance, and connected pathways represent genes within it. C: Correlation network of genes in yellow module, L. crocea and P. plecoglossicida genes with edge weights greater than 0.35. Dot size represents KME values. L. crocea genes are shown in purple, P. plecoglossicida genes are shown in green.

To determine the importance of the modules for our research, we performed enrichment analysis on each module. The yellow module was enriched in many immune-related pathways, as shown in Figure 8B. Therefore, the yellow module was believed to be most important for the infection process. The pathogenic genes in this module were enriched in two meaningful pathways, i.e., flagellar assembly and bacterial chemotaxis.

To understand host-pathogen interactions, we selected several key genes for analysis. We calculated the KME values of genes in the yellow module and selected genes with KME>0.8 and weight>0.35 for visualization. The interaction diagram is shown inFigure 8C. Only five bacterial genes qualified: i.e., dadA, aldH, dapA, L321_RS03170, and L321_RS20070.

DISCUSSION

The low load of pathogens is one of the main limiting factors in the application of dual RNA-seq (Luo et al., 2019b). The spleen is not only an important immune organ for fish (Chen et al., 2019) but also a major target organ for P. plecoglossicida infection (Zhang et al., 2018a). The pathogen load of P. plecoglossicida in the spleen is two or three orders of magnitude higher than that of other organs (Luo et al., 2020). The P. plecoglossicida-infected spleen was chosen for dual RNA-seq, not only because of its high pathogen load, thus meeting the requirements of dual RNA-seq, but also because it can reflect interactions between host and pathogens. In addition, the infected spleen has been used in previous studies on interactions between P. plecoglossicida and host (Sun et al., 2019b; Tang et al., 2019b).

Both the heat map of the major DEMs and PCA analysis showed that the transcriptome changes during infection were enormous. The number of DEMs plot shows the fluctuations and consistency of changes in L. crocea and P. plecoglossicida during infection. The changes in each period of infection were different, as shown in Figures 3, 6. Immune-related GO terms were enriched in the host after 3 dpi, whereas GO functional enrichment results for pathogens differed every day. As infection time increased, the immune system became more sensitive and the immune response was more comprehensive.

As life activity is a dynamic process, time series experiments can effectively characterize the expression of genes. Different from traditional expression difference analysis, it is more crucial to screen expression patterns of interest based on cluster analysis of gene expression patterns, that is, variation in expression levels and subsequent functional enrichment analysis on genes in the given module.

Host innate immunity uses a variety of recognition receptors, such as Toll-like receptors, to identify bacteria that attack the host (Gulati et al., 2018). Transcriptome analysis of Oreochromis niloticus during Streptococcus agalactiae infection revealed that the Toll-like receptor-mediated pathway helps the immune response and protects the host against pathogens (Ken et al., 2017). Here, results showed that the expression of many genes in the Toll-like receptor signaling pathway was significantly up-regulated (Figure 3C) at the later stage of infection (4 dpi). Cytokines are key regulators of host defense in innate immunity and adaptive inflammation (Akdis et al., 2016). Gene expression within the “cytokine-cytokine receptor interaction” changed significantly at 3 dpi and 4 dpi (Figure 3C). Gene expression in module profile 43-LC, which was enriched in the “cytokine-cytokine receptor interaction” pathway based on STEM cluster analysis, was higher than that in the CG overall during infection (Figure 4). These results indicate that the immune response of L. crocea was enhanced during P. plecoglossicida infection.

A host’s metabolic regulation during infection is also very important and has a significant impact on the immune response (Nhan et al., 2019). The immune system can resist pathogens and maintain tissue homeostasis for the life of the organism but requires a large investment in bioenergy (Ganeshan & Chawla, 2014). Thus, to overcome the stresses caused by pathogens, infected fish need to consume large amounts of energy (Yin et al., 2014, 2015). Amino acid metabolism is essential for consumption of energy, synthesis of detoxified protein, and safe operation of innate immune responses in fish (Lu et al., 2017; Zhang et al., 2019b). As seen from STEM analysis, gene expression in profile 18-LC decreased slightly on the first day, but overall showed an upward trend. Based on KEGG enrichment, the genes in profile 18-LC were closely related to the tryptophan metabolism pathway. We speculate that on the first day, after the fish were switched to a new environment, they reduced their consumption to adapt to the environment, accompanied by increasingly serious infection, requiring amino acid metabolism to provide required elements such as energy and protein. Tryptophan and its metabolites play important roles in regulating immune mechanisms in marine animals, especially humoral immunity (Curti et al., 2009; Mondanelli et al., 2019). Lu et al. (2017) claimed that Vibrio parahaemolyticus infection can cause changes in tryptophan metabolism in female Haliotis diversicolor, with similar phenomena also observed in zebrafish (Ji et al., 2019). Deeper infection with Aeromonas salmonicida is also known to reduce metabolites (including choline, glycerophosphocholine (GPC), betaine, and glycine) in the tryptophan metabolism pathway of Atlantic salmon (Liu et al., 2016).

Respiratory bursts are a basic strategy for fish species to fight against pathogens and are often accompanied by substantial oxygen consumption (Biller & Takahashi, 2018; Lu & Chen, 2019; Yilmaz, 2019). During fish immunological responses to pathogens, respiratory bursts play an important biochemical role that facilitates successful phagocytosis and eventual pathogen destruction (Lulijwa et al., 2019). In the current study, however, profile 0-LC, which showed down-regulated expression during infection, exhibited functional enrichment associated with oxygen metabolism (GO: 0005344, GO: 0019825, GO: 0020037, GO: 0005506) (Figure 4, 5); in other words, L. crocea was in a state of reduced breathing. Thus, P. plecoglossicida infection inhibited L. crocea respiration, which is consistent with the experimental phenomenon we observed. We speculate that this effect may be one of the main causes of death from visceral white spot disease.

For pathogens invading a host, flagella play a series of key roles, including reaching the optimal host site, colonization or invasion, maintenance at infection site, and post-infection dispersal (Chaban et al., 2015). In this study, the invasion of pathogens was closely related to motility. The functional enrichment results of profile 43-PP showed that its genes were closely related to flagella (Figure 7). During infection, the genes in this module showed an upward trend until the fourth day, after which gene expression decreased slightly due to increased host death. Comprehensive GO enrichment and STEM results suggest that flagella are critical for pathogens in host infection. In the late stage of infection, pathogens are widely distributed in the host; therefore, the role of flagella may be slightly decreased, resulting in a decrease in the expression of related genes. Flagellar motility is an essential process in many phases of a pathogen’s life cycle, and virulence and motility are often intimately linked by complex regulatory networks (Josenhans & Suerbaum, 2002). To date, flagella are known to be involved in the virulence of many aquatic pathogenic bacteria, such as Yersinia ruckeri (Jozwick et al., 2019), Edwardsiella tarda (Morimoto et al., 2019), Aliivibrio salmonicida (Nørstebø et al., 2017), and Vibrio anguillarum (Liu et al., 2014). In addition, the silencing of flagellum genes fliA and flgM significantly reduces the pathogenicity of P. plecoglossicida (Sun et al., 2019a, 2019b).

Dual RNA-seq analysis evidenced the modulation of P. plecoglossicida and L. crocea during infection. The WGCNA results showed that the relationship between host and pathogen was complex, and five pathogenic genes, i.e., aldH, dapA, dadA, L321_RS03170, and L321_RS20070, showed high correlation with the host immune response. These five genes are all related to the metabolic processes of bacteria (Figure 8C). The product dihydrodipicolinate synthase (DHDPS) encoded by the dapA gene is involved in the diaminopimelate (DAP) pathway, which yields essential metabolic building products, namely, meso-DAP and lysine (Dante et al., 1999; Motoyama et al., 2001). The essential amino acid l-lysine is an important precursor for the synthesis of the peptidoglycan cell wall, housekeeping proteins, and virulence factors of bacteria (Peverelli et al., 2016). DHDPS is the product of an essential bacterial gene and is a promising target for antibiotic development (Griffin et al., 2008). The dapA gene is critical in many bacterial species, including Salmonella typhimurium (Becker et al., 2006), Bacillus subtilis (Kobayashi et al., 2003), Mycobacterium tuberculosis (Shrivastava et al., 2016), and Pseudomonas aeruginosa (Impey et al., 2020). The aldH gene is associated with the aldehyde dehydrogenase (ALDH) superfamily, which is widely found in all kingdoms from archaea to mammals, including humans, and participates in a broad variety of metabolic pathways (Brocker et al., 2013; Jackson et al., 2015; Marchitti et al., 2008). The Pseudomonas genus contains many genes that encode ALDHs (Riveros-Rosas et al., 2019). Annotation of the L321_RS20070 gene has shown that it is also closely related to aldehyde dehydrogenase. The enzyme encoded by dadA is also essential for many bacteria to decompose energy substances (Bardaweel et al., 2011; Li & Lu, 2016).

CONCLUSIONS

In this study, we successfully used tissue dual RNA-seq to simultaneously analyze the dynamics of gene expression changes of both host and pathogen in the context of eukaryotic and prokaryotic biological systems, yielding high-resolution transcriptome data. For L. crocea, both metabolic processes and the use of oxygen had strong impact during infection. WGCNA revealed the relationship between the metabolic processes and the entire immune process. For P. plecoglossicida, motility was the most important life activity during infection, and pathogen invasion exhibited strong dependence on flagella. This study represents the first use of a complex network inference model to infer molecular interspecies interactions between a bacterial pathogen and host (cells and tissues) and successfully predict key disease-causing genes. Application of the tools in more infection systems has dramatically expanded our potential ability to uncover complex pathogenic mechanisms at various periods of in vivo infection.

For future studies: (i) In all previous studies plus the present research, time-resolved dual RNA-seq of tissue studies have provided a first step for application of dual RNA-seq to in vivo models of infection. In the four cases, infected tissues were homogenized prior to RNA extraction and sequencing. This finding means that the host fraction is composed of more than a single cell type. Thus, cell-type-specific tissue dual RNA-seq should be used in the future to increase the resolution of the experiment. (ii) The low bacterial burden in infected tissues and fluctuation of burden during different infectious periods is a limitation for time-resolved dual RNA-seq of tissue. Ongoing progress in sample processing and sequencing technologies provides a high-throughput time point (no limit of infection period) for dual RNA-seq of tissues and provides a more complete picture for understanding the mechanisms underlying bacterial in vivo infection. (iii) The feasibility of time-resolved dual RNA-seq of tissue will pose significant challenges for the bioinformatics analysis of resulting datasets, and future efforts will be needed to exploit their full potential.

Dual RNA-seq was used to identify simultaneous host and pathogen gene expression and their interactions associated with infection time in L. crocea. While gene expression is only one of the many biological processes involved, our findings add to a comprehensive understanding of bacterial diseases in fish and strongly suggest that neither host nor pathogen should be studied in isolation when possible. Many associations between gene expression and infection duration were identified, which should help to elucidate the pathogenesis of visceral white spot disease.

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. and L.X.H. supervised the analyses. G.X., L.M,Z. and W.Q.Z. performed fish experiment. G.X. cultivated the bacteria. Y.T. extracted DNA and RNA. Y.T., L.X.H., Y.X.Q., B.W. and N.L. performed bioinformatics analysis. Y.T. wrote the manuscript with the other authors’ input. Q.P.Y., L.X.H. and Y.Q.S. revised the manuscript. All authors read and approved the final version of the manuscript.