Dual Transcriptomic Analyses Unveil Host–Pathogen Interactions Between Salmonella enterica Serovar Enteritidis and Laying Ducks (Anas platyrhynchos)

Abstract

Salmonella enteritidis (SE) is a pathogen that can readily infect ovarian tissues and colonize the granulosa cell layer such that it can be transmitted via eggs from infected poultry to humans in whom it can cause food poisoning. Ducks are an important egg-laying species that are susceptible to SE infection, yet the host–pathogen interactions between SE and ducks have not been thoroughly studied to date. Herein, we performed dual RNA-sequencing analyses of these two organisms in a time-resolved infection model of duck granulosa cells (dGCs) by SE. In total, 10,510 genes were significantly differentially expressed in host dGCs, and 265 genes were differentially expressed in SE over the course of infection. These differentially expressed genes (DEGs) of dGCs were enriched in the cytokine–cytokine receptor interaction pathway via KEGG analyses, and the DEGs in SE were enriched in the two-component system, bacterial secretion system, and metabolism of pathogen factors pathways as determined. A subsequent weighted gene co-expression network analysis revealed that the cytokine–cytokine receptor interaction pathway is mostly enriched at 6 h post-infection (hpi). Moreover, a number of pathogenic factors identified in the pathogen–host interaction database (PHI-base) are upregulated in SE, including genes encoding the pathogenicity island/component, type III secretion, and regulators of systemic infection. Furthermore, an intracellular network associated with the regulation of SE infection in ducks was constructed, and 16 cytokine response-related dGCs DEGs (including IL15, CD40, and CCR7) and 17 pathogenesis-related factors (including sseL, ompR, and fliC) were identified, respectively. Overall, these results not only offer new insights into the mechanisms underlying host–pathogen interactions between SE and ducks, but they may also aid in the selection of potential targets for antimicrobial drug development.

Introduction

Salmonellosis is an important foodborne illness that primarily occurs as a consequence of consuming eggs contaminated with Salmonella enteritidis (SE) (Gantois et al., 2009; Pijnacker et al., 2019) with more than 75% of SE outbreaks being linked to eggs or products derived therefrom (Braden, 2006). For example, an egg-linked international outbreak of SE infection has been occurring in the European Union for many years, posing a significant risk to consumer health and socioeconomic outcomes (Pijnacker et al., 2019). Although efforts to control SE outbreaks by vaccination and related interventional strategies to date have been made, these approaches have not completely prevented SE prevalence amid poultry flocks. Therefore, better strategies need to be further explored to control SE infection.

During SE infection, both the SE and the host mobilize all available resources to win the life-and-death struggle (Aprianto et al., 2016). One the one hand, SE encodes a number of pathogenic factors, such as PAMPs, including components of the bacterial cell wall–peptidoglycans, lipopolysaccharide, flagellin, and secretion system effectors to attack host innate immune (LaRock et al., 2015). Moreover, SE manipulates the cellular mechanisms of host organisms via pathogen–host interactions (PHIs) in order to take advantage of the capabilities of host cells, leading to infections; subsequently, transovarial transmission to offspring results in the spread of SE through the production of contaminated eggs and fecal matter (Gantois et al., 2009; Jennings et al., 2017). On the other hand, the host response to infection depends on innate immunity in which intrinsic mechanisms are responsible for recognizing and responding to SE challenge. The NFKBIA, IL1B, IL8, and CCL4 genes were consistently induced after endotoxin treatment in chicken (Ciraci et al., 2010). Besides this, the AVD, EXFABP, IRG1, AH221, TRAP6, SAA, and all immunoglobulin genes also played important roles in the course of SE infection in chickens (Matulova et al., 2012). These findings only revealed the mechanism of Salmonella in the process of infection from host or pathogen. However, pathogen infection is a complex process involving the interaction between pathogen and host. Understanding of the host–SE interactions may reveal the underlying mechanism of SE infection, which will lay the theoretical foundation for better strategies to control the SE pathogen.

Ducks, especially Shaoxing and Pekin ducks, are an economically important poultry species, which are susceptible to SE infection (Yan et al., 2008; Cha et al., 2013; Zhang et al., 2019a). Notably, SE has recently grown to be the most common Salmonella serotype to be isolated from laying ducks throughout the world, and its prevalence in countries such as China and the United Kingdom have led to major economic losses while threatening human health and contaminating water (Martelli et al., 2016; Yang et al., 2019; Wang et al., 2020). SE is known to readily infect ovarian tissues and to colonize the granulosa cell layer in poultry (Thiagarajan et al., 1994). However, the PHIs between SE and duck granulosa cells (dGCs) have not been effectively defined. As such, we herein generated a time-resolved model of the SE infection of dGCs to simultaneously assess dynamic changes in gene expression in both duck and SE cells via a dual RNA-seq approach, enabling us to characterize the interactions between these two species and to identify key duck immune genes and SE virulence genes related to these responses.

Materials and Methods

Ethics Statement

The Institutional Animal Care and Use Committee of Yangzhou University (Jiangsu, China) approved all animal studies conducted herein, which were consistent with experimental animal use guidelines and regulations. All ducks were housed in a standard facility that met the requirements outlined in the Laboratory Animal Requirements of Environment and Housing Facilities (GB 14925-2001) publication, and all protocols were consistent with the Jiangsu Administration Rules for Laboratory Animal Use.

Experimental Model Animals and Bacterial Strains

Healthy Shaoxing ducks (Anas platyrhynchos, Chinese native breed, 26 weeks old) that were free of Salmonella were obtained from the National Waterfowl Conservation Farm (Taizhou, China). SE (No. MY1, phage type 4) was isolated from ducks and maintained by the Key Laboratory of Animal Disease and Human Health of Sichuan Province (Deng et al., 2008). The National Center for Medical Culture Collection (Beijing, China) established the phage type and serotype of SE isolates.

dGCs Isolation and SE Infection

Duck granulosa cells were isolated and cultured as in prior studies (Gilbert et al., 1977; Zhang et al., 2019b). Briefly, 10–15 adult prehierarchical follicles were collected from egg-laying ducks under sterile conditions. Yolks and vitelline membranes were removed by rinsing these follicles with Ca²⁺- and Mg²⁺-free PBS, after which tissues were minced into 1–2 mm³ pieces and digested with type II collagenase (1 mg/mL; Sigma, St. Louis, MO, United States) for 5 min at 37°C. Samples were then passed through a 200-μm nylon filter, centrifuged twice for 5 min at 67 × g, washed with M199 media to remove the remaining collagenase and cell debris, followed by resuspension in 3 mL of 50% Percoll, centrifugation at 421 × g for 15 min, and cell layer isolation. The dGCs were then suspended in M199 media (5% fetal calf serum, 2 mmol/L L-glutamine, 5 μg/mL transferrin, 10 μg/mL insulin, 1.75 mM HEPES [4-(2-hydroxyerhyl) piperazine-1-erhanesulfonic acid]) and counted with 0.1% trypan blue being added to assess viability. Only samples with >90% cell survival were used for subsequent experiments. Cells were allowed to adhere in tissue culture flasks for 24 h prior to experimental utilization, after which purity was initially determined by indirect immunofluorescence assay (IFA) to detect follicle granulosa cell-specific receptor, follicle-stimulating hormone receptor (FSHR) expression using anti-FSHR antibodies (1:500, Proteintech, Rosemont, IL, United States, L594-22665).

Invasion Assays

Invasion assays were conducted as in prior reports (Jouve et al., 1997; Lane and Mobley, 2007). Briefly, bacteria were cultured in LB broth to an OD₆₀₀ of 2.0 (at mid logarithmic phase) at 37°C. dGCs were added to 96-well plates (1 × 10⁵/well) for 24 h, after which they were washed thrice with PBS (pH 7.2) prior to the addition to each well of 100 μL of SE in Dulbecco’s modified Eagle medium (DMEM) at a multiplicity of infection (MOI) of 10 with DMEM being added to control wells. Plates were then cultured 37°C for 1 h, after which PBS supplemented with 50 μg/mL gentamicin was added to kill any non-invasive bacteria. Cells were then incubated for an additional 1 h period and then lysed by 0.5% Triton X-100 (Solarbio, Beijing, China) for 30 min before vigorous pipetting. Samples were removed, diluted, and plated on LB agar to determine the number of colony forming units (CFU) per monolayer. Samples were assessed in triplicate.

Salmonella enteritidis infection was confirmed via an IFA approach as detailed previously (Zhang et al., 2019c). Briefly, medium was removed, and cells were washed three times with 0.01 M PBS, fixed with 4% paraformaldehyde, permeabilized with 0.1% Triton X-100, and blocked for 1 h with 10% FBS in PBS. After a 2-h incubation with an anti-Salmonella antibody (1:2000, Abcam, Cambridge, MA, United States, ab69253), cells were washed thrice with PBS, stained with DAPI (0.2 mg/mL) for 15 min at 37°C, and imaged via fluorescence microscopy (Leica, Wetzlar, Germany).

RNA-Sequencing

TRIzol (Invitrogen, Carlsbad, CA, United States) was used based on provided directions to extract RNA from 0, 3, 6, and 9 hpi groups (n = 3, respectively) with 1% agarose gel electrophoresis being used to assess RNA quality and contamination. A NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, United States) was also used to assess RNA purity, and a Qubit RNA Assay Kit and a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States) were used to measure RNA concentrations. An RNA Nano 6000 Assay Kit and a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States) were then used to assess RNA integrity. Next, a strand-specific library was constructed followed by the previous method (Parkhomchuk et al., 2009; Borodina et al., 2011). Briefly, oligo (dT) beads were used to enrich for eukaryotic RNA in each sample first, then prokaryotic mRNA was enriched via the removal of rRNA with a Ribo-Zero Magnetic Kit (Epicentre) and Quick-RNA^TM Fungal/Bacterial Microprep (Zymo Research). After enrichment, fragmentation buffer was used to fragment mRNA into short segments, after which first-strand DNA synthesis was conducted using random hexamer primers. Second-strand cDNA synthesis was then performed with DNA polymerase I, RNase H, dNTPs, and an appropriate buffer. A QiaQuick PCR extraction kit was then used for cDNA isolation, followed by end repair, polyadenylation, and Illumina sequencing adapter ligation. Ligation product size selection was then performed via agarose gel electrophoresis, followed by PCR amplification and sequencing runs of paired-end 150 bp with an Illumina HiSeq 2500 instrument at Major BioTech Co., Ltd. (Shanghai, China).

Raw reads were trimmed and subjected to quality control using Sickle¹ and SeqPrep² with default settings. After cleaning, data were mapped to the duck [NCBI: IASCAAS_PekingDuck_PBH1.5 (GCF_003850225.1)] and SE (NC_011294.1) reference genomes using TopHat2 and Bowtie2, respectively (Trapnell et al., 2010; Langmead and Salzberg, 2012). The “MarkDuplicate” tools of Picard³ were used removed the reads that mapped to both organisms. Relationships among samples were assessed with the gmodels R package⁴.

Differentially Expressed Gene (DEG) Analyses

The gene abundances were quantified by software RSEM, and the gene expression level was normalized by using the fragments per kilobase of transcript per million mapped reads method (Li and Dewey, 2011). Differentially expressed genes (DEGs) among samples were identified with the edgeR package (Robinson et al., 2010). DEGs were defined as those genes with a | log2FC| ≥ 2 and a false discovery rate (FDR) < 0.05. After identification, Blast2GO (Conesa and Götz, 2008) was utilized for the GO annotation of these genes (Conesa et al., 2005), and the KEGG program was utilized for KEGG pathway enrichment analyses (Kanehisa and Goto, 2000). The weighted gene co-expression network analysis (WGCNA) R package was used to conduct a WGCNA of gene expression patterns across samples (Langfelder and Horvath, 2008). Following pair-wise correlation matrix calculation, an adjacency matrix was calculated by raising the correlation matrix based on a power value (β) that had been selected based upon scale-free topology criteria (scale-free topology criteria R² ≥ 0.85). Topological overlap was then assessed. A cluster tree branch gene co-expression module was constructed based upon average linkage hierarchical clustering. Module eigengenes were statistically classified (dynamicTreeCut algorithm, parameters deepSplit = 2, minModulesize = 30) (Langfelder et al., 2008), after which very similar modules were merged, leading to the identification of key host genes associated with responses to SE infection at different time points. Clusterprofiler and DAVID were used for KEGG and GO enrichment analysis, respectively (Huang et al., 2009; Yu et al., 2012).

Interaction Network Construction and Prediction of the Regulation of Pathogenicity-Associated Genes

Genes involved in PHIs were identified by using BLASTp to compare potentially relevant genes against the PHI-base version 4.8 database (September 2019) (Urban et al., 2020). Network-based inferences necessitating a gene expression profile is centralized and standardized using the min-max normalization method (Zhu et al., 2019). Gene expression pattern analysis is used to cluster genes of similar expression patterns for multiple samples. To examine the expression pattern of DEGs, the expression data of each sample (in the order of treatment) were normalized to 0, log2 (v1/v0), log2 (v2/v0), and then clustered by Short Time-series Expression Miner software (STEM) (Ernst and Bar-Joseph, 2006). Host dGCs DEGs (euk-DEGs) and prokaryotic SE-associated DEGs (pro-DEGs) were then selected to predict associations between key virulence genes and host genes via an expression correlation analysis. Interaction networks of euk-DEGs and corresponding pro-DEGs were constructed using Cytoscape (Shannon et al., 2003). To improve network visibility, the edges (connections) were filtered to show only those with a degree cutoff = 2, a node score cutoff = 0.2, k-core = 2, and max depth = 100. The two significant networks were then calculated with the MCODE plug-in (Bader and Hogue, 2003).

qPCR-Based Result Validation

In order to validate RNA-seq results via qPCR, 12 duck genes (IL15, CCR7, IL18R1, CD40, IL23R, IL20RA, TNFRSF8, TNFRSF13B, TNFRSF4, CSF3R, IL2RG, TNFRSF18) and 17 SE virulence-associated genes (sseL, sseE, invE, ssaM, ssaN, ssaP, spiC, slyA, spaR, aroA, sipB, prgH, hilD, ompR, fruR, fliC, sipA) were selected for analysis with GAPDH and gyrA being used for the respective normalization of duck and SE transcript expression levels. Primers were designed with the Primer 5 software (Premier, Canada) based upon sequencing data (Table 1), and were synthesized by Sangon Biotech (Shanghai, China). TRIzol (Takara, Dalian, China) was used based on provided directions to extract RNA from prepared samples. All qPCR reactions were conducted with an ABI 7500 instrument (Applied Biosystems, Foster City, CA, United States) with individual reactions being composed of 10 μL 2 × SYBR Premix Ex Taq II (Takara), 0.4 μL of each primer (10 μM), 0.4 μL 50 × ROX Reference Dye II (Takara), 2 μL cDNA, and ddH₂O to a final volume of 20 μL. Thermocycler settings were 95°C for 30 s, 40 cycles of 95°C for 5 s, 60°C for 34 s. The thermal program for RT-qPCR was 98°C for 30 s, followed by 40 cycles of 98°C for 10 s and 60°C for 30 s. All samples were conducted with three technical replications. The 2^–ΔΔCT method (Kenneth and Thomas, 2002) was used to assess relative gene expression levels.

Statistical Analyses

Data are means ± standard deviation (SD), and experiments were repeated at least three times. Data were analyzed with SPSS 26.0 (Chicago, IL, United States) via one-way ANOVAs with Dunnett’s test, with p < 0.05 as the significance threshold.

Data Access

RNA-seq data from the present study have been deposited in the Short Read Archive of the National Center for Biotechnology Information (NCBI) under the bio-project number PRJNA721551.

Results

Evaluation of dGC Invasion by SE

After isolation, dGCs exhibited a rounded or ovoid shape and were fully adherent to tissue culture flasks within 24 h of culture at which point they were spread out in a pebble-like manner. These cells reached 100% confluence within 3–5 days, at which time granular material was evident in the cytoplasm. DAPI was used to stain dGCs nuclei, and the cytoplasm of the FSHRs was stimulated to emit red fluorescence for detection by IFA (Figure 1A). After treatment with bacterial suspensions and gentamicin-containing growth media, changes in cell morphology and numbers of invasive SE were assessed at different time points post-infection. At 3 hpi, cells appeared swollen and rounded or irregular in shape with some evidence of vacuolation. At 6 hpi, infected cells were increasingly non-adherent with some blue or black particles, cellular deformation, and clumping in strip-like aggregates. At 9 hpi, some cells had broken down and dissolved with increases in the presence of dark blue particles within dGCs. The media also appeared increasingly flocculated at this time point due to cell exfoliation. IFA staining of the SE in these samples revealed that the degree of bacterial invasion increased significantly over time such that relatively few bacteria were detectable at 3 hpi, whereas these bacterial numbers had expanded significantly by 6 hpi and had invaded most cells by 9 hpi (Figure 1B). The intracellular CFU counts of SE at each time point are listed in Supplementary Table 1.

FIGURE 1

Profiling of Gene Expression Patterns Associated With HPIs

To examine dGCs and SE gene expression profiles over the course of this infection process, we conducted a dual RNA-seq time-course analysis of both dGCs and SE cells at 0, 3, 6, and 9 hpi (Supplementary Table 2). This RNA-seq approach yielded high-quality yields with an appropriately balanced base distribution and a reasonable N% (>10%) consistent with these data being of sufficient quality to permit downstream analyses (Supplementary Table 3). Average base error rate distributions in these sequencing reads were <0.1% and were, thus, considered acceptable. Principal component analysis clustering was then conducted to confirm sample distribution and biological reproducibility (Figures 2A,B). The gmodels R package was used to evaluate the relationships among samples (Supplementary Figure 1). On average, this analysis yielded 78.6 million reads per library with 76.2% of these reads aligning to the duck reference genome and 23.8% aligning to the SE P125109 reference genome (Supplementary Table 4 and Figure 2C). The edgeR package (Robinson et al., 2010) was used to assess gene expression profiles and to identify significant DEGs that met the defined criteria (FDR < 0.05 and | log2FC| > 1) at different time points in both dGCs and SE cells (Figure 2D). For further details regarding DEGs and virulence-associated genes between the 0, 3, 6, and 9 hpi time points and associated biological functions, see Supplementary Tables 5, 6.

FIGURE 2

Enrichment Analysis of Pathways and GO Terms for DEGs in SE and dGCs

GO analyses of identified DEGs in dGCs and SE were conducted (Supplementary Tables 7, 8), revealing these genes to be associated with GO terms mainly enriched in host immune function, including immune system process (GO: 0002376), regulation of cytokine production (GO: 0001817), positive regulation of biological process (GO: 0001819), and regulation of cellular metabolic process (GO: 0031323). Bacterial DEGs were additionally associated with bacterial virulence–related protein secretion, such as the metabolic process (GO: 0008152), peptide biosynthetic process (GO: 0043043), and regulation of gene expression (GO: 0010468) GO terms (Supplementary Figure 2). In addition, the top 20 enriched KEGG pathways associated with genes that were differentially expressed during different stages of the SE infection process were assessed (Figure 3). These pathways included the cytokine–cytokine receptor interaction, Toll-like receptor signaling, JAK-STAT signaling, ECM-receptor interaction, and p53 signaling pathways, which play key roles in immune, signal transduction, and transcription-related processes. Statistical analyses suggest that the cytokine–cytokine receptor interaction pathway in dGCs was significantly altered during SE infection as were the two-component system and bacterial secretion system pathways in SE (Table 2). For further details regarding the DEGs associated with the cytokine–cytokine receptor interaction, two-component system, and bacterial secretion system pathways in this experimental system, see Table 3.

FIGURE 3

Identification of Key Genes Associated With dGC Responses to Infection and SE Pathogenicity

A WGCNA analysis was conducted to assess relationships between particular genes during a given stage of the infection process, grouping similarly expressed genes into modules via average linkage clustering (based on the weighted correlation coefficients of genes, genes are classified according to their expression patterns, and genes with similar patterns are grouped into one module). Herein, a power value of β = 9 (scale-free R² = 0.85) was utilized for oft-thresholding to ensure a scale-free network (Figure 4A), leading to the identification of 12 modules (Figure 4B). Module connectivity was assessed based upon Pearson’s correlation coefficient values (cor.geneModuleMembership > 0.8) with clinical trait relationships similarly being measured based upon absolute Pearson’s correlation coefficient values (cor.geneTraitSignificance > 0.4). As the identification of modules significantly associated with clinical features was considered to be of biological importance, module–feature relationships were assessed at different pathological stages, revealing that the most closely associated modules at the 0, 3, 6, and 9 hpi time points were the blue, dark red, black, and light cyan modules, respectively (Figures 4C,D and Supplementary Table 9). These four modules were, thus, selected as modules of interest with clinical features worthy of examination in subsequent analyses. Genes in these significant modules at the 0, 3, 6, and 9 hpi time points were associated with metabolic pathways, the MAPK and JAK-STAT signaling pathways, the cytokine–cytokine receptor interaction pathway, and the cell cycle and p53 signaling pathways, respectively (Table 2 and Supplementary Figure 3). The key genes in the cytokine–cytokine receptor interaction pathway within the black modules are additionally listed in Table 4. To further understand SE pathogen-related factor expression over the course of the infection process, identified genes were subjected to a BLAST analysis and annotation of high-abundance genes (expression > 10) using the PHI database with identified pathogenicity-associated factors being listed in Table 4.

FIGURE 4

Co-expression Analysis of DEGs in SE and dGCs

Fuzzy c-means clustering (Bezdek, 1992) and expression correlation analyses were next used to evaluate the expression matrices of DEGs from dGCs and SE cells (Supplementary Figures 4,5). This led to the identification of 16 host genes associated with the cytokine–cytokine receptor interaction pathway that exhibited high fold-change expression patterns during the infection process as well as 17 pathogenicity-associated factors that were primarily associated with functions including regulators of systemic infection, pathogenicity island/component, and type III secretion system components (<0.05, P_x,y > 0.6) (Figure 5). In Figure 5, host and pathogen DEGs are, respectively, represented by circles and inverted triangles. The expression levels of most host cytokines were positively correlated with the expression of bacterial virulence factors, particularly the type III secretion-related proteins; bacterial virulence factors ompR, fliC, and ssaU were the negative regulators to the expression of host response genes; and sseL, ttrB, and ompR interact throughout the higher cooperative diversity in the regulatory network.

FIGURE 5

qPCR-Based Validation of Hub Genes and Pathogenicity-Associated Factors

Last, a qPCR approach was used to validate the DEGs identified in the above RNA-seq analysis. In total, 12 hub genes (IL15, CCR7, IL18R1, CD40, IL23R, IL20RA, TNFRSF8, TNFRSF13B, TNFRSF4, CSF3R, IL2RG, and TNFRSF18) and 17 pathogenicity-associated factors (sseL, sseE, invE, ssaM, ssaN, ssaP, spiC, slyA, spaR, aroA, sipB, prgH, hilD, ompR, fruR, fliC, and sipA) were selected for validation. Subsequent qPCR results appeared similar to the outcomes of the dual RNA-seq analysis, thus confirming the validity of these RNA-Seq results (Figure 6).

FIGURE 6

Discussion

Salmonella enteritidis is a zoonotic pathogen that can readily colonized the granulosa cell layer of ovarian tissues (Gantois et al., 2009). Clarifying HPIs will help reveal the mechanism for regulating the virulence of pathogenic bacteria. In this study, the 10,510 DEGs in host dGCs and 265 in SE were identified during SE infection, and 16 cytokine response-related dGCs DEGs (including IL15, CD40, and CCR7) and 17 pathogenesis-related factors (including sseL, ompR, and fliC) were identified. These results were similar to those of chicken granulosa cells or duck ovary tissues challenged by SE (Tsai et al., 2010; Babu et al., 2016; Zhang et al., 2019a). Other genes, including IAP1, CD28, TGFβ2-4, Gal11-13, TRAIL, PSAP, CASP1, AVD, EXFABP, IRG1, and AH221, et al., played an important role of anti-SE immune regulation in cecum and spleen (Calenge and Beaumont, 2012; Matulova et al., 2012, 2013). Together, these data suggest that the immune response of host immune genes were different in different tissues and organs during SE infection. Salmonella, as a pathogen of vertical transmission through the ovary, should establish an ovary-specific immune response during the infection process.

Furthermore, the cytokine–cytokine receptor interaction pathway, two-component system, and bacterial secretion system were enriched to involve HPIs. Results of the present study were not in agreement with the findings of Huang et al. (2016) and Wang et al. (2019), who demonstrate that the IgA production signaling pathway and TLR4-FOS/JUN pathway contributed to the protection of chicken from Salmonella invasion, respectively. The reason for these different pathways might be due to sequencing methods. In the present study, the dual transcriptomic analyses were performed to unveil HPIs. Westermann and Vogel (2018) prove that dual RNA sequencing is an effective strategy to reveal the pathogenic mechanisms and the host immune response although the RNA sequencing could only reflect the immune response of the pathogen or host. In this sense, the cytokine–cytokine receptor interaction pathway, two-component system, and bacterial secretion system are more likely to be involved in the immune response during Salmonella infection.

Evaluating host–SE interactions not only helps uncover the pathogenesis of SE infections, but also has the potential to develop to antimicrobial drug targeting. For example, sMtb-RECON was established by a combination model of bacteria Mycobacterium tuberculosis (Mtb) and human macrophage-like cell line THP-1 metabolic processes used as the elucidation of metabolic drug responses in a manner that has the potential to reduce antibiotic abuse (Rienksma et al., 2019). In another study, a dual RNA-sequencing time-course approach was used to monitor transcriptomic responses in both Salmonella typhimurium and infected HeLa cells, revealing a link between the S. typhimurium PinT virulence gene and the induction of immune signaling responses in HeLa cells via JAK-STAT signaling pathways (Westermann et al., 2016). In the present study, we find that BPSS1504 was involved in SE infection, which served as a Type VI secretion system (T6SS) component capable of influencing hemolysin-coregulated protein Hcp1 secretion and T6SS apparatus integrity (Hopf et al., 2014). The interaction between BPSS1504 and T6SS might be a potential target to develop an antimicrobial drug for SE clearance.

Conclusion

In conclusion, the 10,510 euk-DEGs and 265 pro-DEGs were screened during duck SE challenge using dual RNA-seq, which enriched in the host cytokine–cytokine receptor interaction pathway and bacterial two-component/secretion system throughout the process, respectively. Also, a number of PHIs between cytokine–cytokine receptors and virulence-associated genes were found in this process based on WGCNA. Furthermore, the intracellular network associated with the regulation of SE infection in ducks was constructed, and 16 cytokine response–related dGC DEGs (including IL15, CD40, and CCR7, et al.) and 17 pathogenesis-related factors (including sseL, ompR, and fliC, et al.) were identified. Also, the interaction between BPSS1504 and T6SS might be involved in immune regulation during SE infection. These data offer novel insights into the complex HPIs that occur upon SE infection of ducks, and unveil novel potential targets for SE infection diagnosis and antimicrobial drug development.

Publisher’s Note

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

References

• 1

  ApriantoR.SlagerJ.HolsappelS.VeeningJ. W. (2016). Time-resolved dual RNA-seq reveals extensive rewiring of lung epithelial and pneumococcal transcriptomes during early infection.Genome Biol.17:198. 10.1186/s13059-016-1054-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  BabuU. S.HarrisonL. M.PatelI. R.RamirezG. A.WilliamsK. M.PereiraM.et al (2016). Differential antibacterial response of chicken granulosa cells to invasion by Salmonella serovars.Poultry Sci.951370–1379. 10.3382/ps/pew050

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  BaderG. D.HogueC. W. (2003). An automated method for finding molecular complexes in large protein interaction networks.BMC Bioinform.4:2. 10.1186/1471-2105-4-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  BezdekJ. C. (1992). “Clustering in banach spaces”, in Fuzzy Approach to Reasoning and Decision-Making, ed.NovákV. (Dordrecht: Springer), 173–184. 10.1007/978-94-011-2753-0_15

  • CrossRef
  • Google Scholar

• 5

  BorodinaT.AdjayeJ.SultanM. (2011). A strand-specific library preparation protocol for RNA sequencing.Methods Enzymol.50079–98. 10.1016/B978-0-12-385118-5.00005-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  BradenC. R. (2006). Salmonella enterica serotype Enteritidis and eggs: a national epidemic in the United States.Clin. Infect. Dis.43512–517. 10.1086/505973

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  CalengeF.BeaumontC. (2012). Toward integrative genomics study of genetic resistance to Salmonella and Campylobacter intestinal colonization in fowl.Front. Genet.3:261. 10.3389/fgene.2012.00261

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  ChaS. Y.KangM.YoonR. H.ParkC. K.MoonO. K.JangH. K. (2013). Prevalence and antimicrobial susceptibility of Salmonella isolates in Pekin ducks from South Korea.Comp. Immuno. Microbiol. Infect. Dis.36473–479. 10.1016/j.cimid.2013.03.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  CiraciC.TuggleC. K.WannemuehlerM. J.NettletonD.LamontS. J. (2010). Unique genome-wide transcriptome profiles of chicken macrophages exposed to Salmonella-derived endotoxin.BMC Genom.11:545. 10.1186/1471-2164-11-545

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  ConesaA.GötzS. (2008). Blast2GO: A comprehensive suite for functional analysis in plant genomics.Int. J. Plant Genom.2008:619832. 10.1155/2008/619832

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  ConesaA.GotzS.Garcia-GomezJ. M.TerolJ.TalonM.RoblesM. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research.Bioinformatics.213674–3676. 10.1093/bioinformatics/bti610

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  DengS. X.ChengA. C.WangM. S.CaoP. (2008). Serovar-specific real-time quantitative detection of Salmonella enteritidis in the gastrointestinal tract of ducks after oral challenge.Avian Dis.5288–93. 10.1637/8102-090107-Reg

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  ErnstJ.Bar-JosephZ. (2006). STEM: a tool for the analysis of short time series gene expression data.BMC Bioinform.7:191. 10.1186/1471-2105-7-191

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  GantoisI.DucatelleR.PasmansF.HaesebrouckF.GastR.HumphreyT. J.et al (2009). Mechanisms of egg contamination by Salmonella Enteritidis.FEMS Microbiol. Rev.33718–738. 10.1111/j.1574-6976.2008.00161.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  GilbertA. B.EvansA. J.PerryM. M.DavidsonM. H. (1977). A method for separating the granulosa cells, the basal lamina and the theca of the preovulatory ovarian follicle of the domestic fowl (Gallus domesticus).J. Reprod. Fertil.50179–181. 10.1530/jrf.0.0500179

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  HopfV.GöhlerA.Eske-PogoddaK.BastA.SteinmetzI.BreitbachK. (2014). BPSS1504, a cluster 1 type VI secretion gene, is involved in intracellular survival and virulence of Burkholderia pseudomallei.Infect. Immun.822006–2015. 10.1128/IAI.01544-14

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  HuangD. W.ShermanB. T.LempickiR. A. (2009). Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists.Nucleic Acids Res.371–13. 10.1093/nar/gkn923

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  HuangH. B.LiuA.WuH.AnsariA. R.WangJ. X.HuangX. Y.et al (2016). Transcriptome analysis indicated that Salmonella lipopolysaccharide-induced thymocyte death and thymic atrophy were related to TLR4-FOS/JUN pathway in chicks.BMC Genom.17:322. 10.1186/s12864-016-2674-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  JenningsE.ThurstonT. L. M.HoldenD. W. (2017). Salmonella SPI-2 type III secretion system effectors: molecular mechanisms and physiological consequences.Cell. Host Microb.22217–231. 10.1016/j.chom.2017.07.009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  JouveM.GarciaM. I.CourcouxP.LabigneA.GounonP.Le BouguénecC. (1997). Adhesion to and invasion of HeLa cells by pathogenic Escherichia coli carrying the afa-3 gene cluster are mediated by the AfaE and AfaD proteins, respectively.Infect. Immun.654082–4089. 10.1128/iai.65.10.4082-4089.1997

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  KanehisaM.GotoS. (2000). KEGG: kyoto encyclopedia of genes and genomes.Nucleic Acids Res.2827–30. 10.1093/nar/28.1.27

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  KennethJ. L.ThomasD. S. (2002). Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method.Methods25402–408. 10.1006/meth.2001.1262

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  LaneM. C.MobleyH. L. T. (2007). Role of P-fimbrial-mediated adherence in pyelonephritis and persistence of uropathogenic Escherichia coli (UPEC) in the mammalian kidney.Kidney Int.7219–25. 10.1038/sj.ki.5002230

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  LangfelderP.HorvathS. (2008). WGCNA: an R package for weighted correlation network analysis.BMC Bioinform.9:559. 10.1186/1471-2105-9-559

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  LangfelderP.ZhangB.HorvathS. (2008). Defining clusters from a hierarchical cluster tree: the dynamic tree cut package for R.Bioinformatics.24719–720. 10.1093/bioinformatics/btm563

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  LangmeadB.SalzbergS. L. (2012). Fast gapped-read alignment with Bowtie 2.Nat. Methods.9357–359. 10.1038/nmeth.1923

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  LaRockD. L.ChaudharyA.MillerS. I. (2015). Salmonellae interactions with host processes.Nat. Rev. Microbiol.13191–205. 10.1038/nrmicro3420

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  LiB.DeweyC. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome.BMC Bioinform.12:323. 10.1186/1471-2105-12-323

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  MartelliF.BirchC.DaviesR. H. (2016). Observations on the distribution and control of Salmonella in commercial duck hatcheries in the UK.Avian Pathol.45261–266. 10.1080/03079457.2016.1146820

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  MatulovaM.RajovaJ.VlasatikovaL.VolfJ.StepanovaH.HavlickovaH.et al (2012). Characterization of chicken spleen transcriptome after infection with Salmonella enterica Serovar Enteritidis.PloS One.7:e48101. 10.1371/journal.pone.0048101

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  MatulovaM.VarmuzovaK.SisakF.HavlickovaH.BabakV.StejskalK.et al (2013). Chicken innate immune response to oral infection with Salmonella enterica serovar Enteritidis.Vet. Res.44:37. 10.1186/1297-9716-44-37

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  ParkhomchukD.BorodinaT.AmstislavskiyV.BanaruM.HallenL.KrobitschS.et al (2009). Transcriptome analysis by strand-specific sequencing of complementary DNA.Nucleic Acids Res.37:e123. 10.1093/nar/gkp596

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  PijnackerR.DallmanT. J.TijsmaA. S. L.HawkinsG.LarkinL.KotilaS. M.et al (2019). An international outbreak of Salmonella enterica serotype Enteritidis linked to eggs from Poland: a microbiological and epidemiological study.Lancet Infect. Dis.19778–786. 10.1016/S1473-3099(19)30047-7

  • CrossRef
  • Google Scholar

• 34

  RienksmaR. A.SchaapP. J.dos SantosV. A. P. M.Suarez-DiezM. (2019). Modeling host-pathogen interaction to elucidate the metabolic drug response of intracellular Mycobacterium tuberculosis.Front. Cell. Infect. Microbiol.9:144. 10.3389/fcimb.2019.00144

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.Bioinformatics.26139–140. 10.1093/bioinformatics/btp616

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD.et al (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks.Genome Res.132498–2504. 10.1101/gr.1239303

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  ThiagarajanD.SaeedA. M.AsemE. K. (1994). Mechanism of transovarian transmission of Salmonella enteritidis in laying hens.Poult. Sci.7389–98. 10.3382/ps.0730089

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  TrapnellC.WilliamsB. A.PerteaG.MortazaviA.KwanG.van BarenM. J.et al (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation.Nat. Biotech.28511–515. 10.1038/nbt.1621

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  TsaiH. J.ChiuC. H.WangC. L.ChouC. H. (2010). A time-course study of gene responses of chicken granulosa cells to Salmonella Enteritidis infection.Vet. Microbiol.144325–333. 10.1016/j.vetmic.2010.01.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  UrbanM.CuzickA.SeagerJ.WoodV.RutherfordK.VenkateshS. Y.et al (2020). PHI-base: the pathogen-host interactions database.Nucleic Acids Res.48D613–D620. 10.1093/nar/gkz904

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  WangF.ZhangJ.ZhuB.WangJ.WangQ.ZhengM. Q.et al (2019). Transcriptome analysis of the Cecal Tonsil of Jingxing yellow chickens revealed the mechanism of differential resistance to Salmonella.Genes10:979. 10.3390/genes10120979

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  WangJ.LiJ. X.LiuF. L.ChengY. Y.SuJ. L. (2020). Characterization of Salmonella enterica Isolates from diseased poultry in Northern China between 2014 and 2018.Pathogens9:95. 10.3390/pathogens9020095

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  WestermannA. J.ForstnerK. U.AmmanF.BarquistL.ChaoY. J.SchulteL. N.et al (2016). Dual RNA-seq unveils noncoding RNA functions in host-pathogen interactions.Nature529496–501. 10.1038/nature16547

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  WestermannA. J.VogelJ. (2018). Host-pathogen transcriptomics by dual RNA-Seq.Methods Mol. Biol.173759–75. 10.1007/978-1-4939-7634-8_4

  • CrossRef
  • Google Scholar

• 45

  YanB.ChengA. C.WangM. S.DengS. X.ZhangZ. H.YinN. C.et al (2008). Application of an indirect immunofluorescent staining method for detection of Salmonella enteritidis in paraffin slices and antigen location in infected duck tissues.World J. Gastroenterol.14776–781. 10.3748/wjg.14.776

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  YangJ.JuZ.YangY.ZhaoX.JiangZ.SunS. (2019). Serotype, antimicrobial susceptibility and genotype profiles of Salmonella isolated from duck farms and a slaughterhouse in Shandong province, China.BMC Microbiol.19:202. 10.1186/s12866-019-1570-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  YuG. C.WangL. G.HanY. Y.HeQ. Y. (2012). clusterProfiler: an R Package for comparing biological themes among gene clusters.Omics16284–287. 10.1089/omi.2011.0118

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  ZhangY.ChenY.GuT.XuQ.ZhuG.ChenG. (2019a). Effects of Salmonella enterica serovar enteritidis infection on egg production and the immune response of the laying duck Anas platyrhynchos.PeerJ.7:e6359. 10.7717/peerj.6359

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  ZhangY.GuT.ChenY.ZhuG.VongsangnakW.XuQ.et al (2019b). Screening and identification of SipC-interacting proteins in Salmonella enteritidis using Gal4 yeast two-hybrid system in duck.PeerJ7:e7663. 10.7717/peerj.7663

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  ZhangY.GuT. T.ChenY.HuangY.DuJ. P.LuL. Z.et al (2019c). Comparative transcriptome analysis reveals PERP upregulated during Salmonella Enteritidis challenge in laying ducks.J. Cell. Physiol.23411330–11347. 10.1002/jcp.27790

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  ZhuX. L.FuB. F.YangY. D.MaY.HaoJ. Y.ChenS. Q.et al (2019). Attention-based recurrent neural network for influenza epidemic prediction.BMC Bioinform.20:575. 10.1186/s12859-019-3131-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar