Skip to main content
  • Research article
  • Open access
  • Published:

Transcriptome analysis reveals modulation of the STAT family in PEDV-infected IPEC-J2 cells

Abstract

Background

Porcine epidemic diarrhea virus (PEDV) is a causative agent of serious viral enteric disease in suckling pigs. Such diseases cause considerable economic losses in the global swine industry. Enhancing our knowledge of PEDV-induced transcriptomic responses in host cells is imperative to understanding the molecular mechanisms involved in the immune response. Here, we analyzed the transcriptomic profile of intestinal porcine epithelial cell line J2 (IPEC-J2) after infection with a classical strain of PEDV to explore the host response.

Results

In total, 854 genes were significantly differentially expressed after PEDV infection, including 716 upregulated and 138 downregulated genes. Functional annotation analysis revealed that the differentially expressed genes were mainly enriched in the influenza A, TNF signaling, inflammatory response, cytokine receptor interaction, and other immune-related pathways. Next, the putative promoter regions of the 854 differentially expressed genes were examined for the presence of transcription factor binding sites using the MEME tool. As a result, 504 sequences (59.02%) were identified as possessing at least one binding site of signal transducer and activator of transcription (STAT), and five STAT transcription factors were significantly induced by PEDV infection. Furthermore, we revealed the regulatory network induced by STAT members in the process of PEDV infection.

Conclusion

Our transcriptomic analysis described the host genetic response to PEDV infection in detail in IPEC-J2 cells, and suggested that STAT transcription factors may serve as key regulators in the response to PEDV infection. These results further our understanding of the pathogenesis of PEDV.

Background

Porcine epidemic diarrhea (PED), caused by the porcine epidemic diarrhea virus (PEDV), is one of the most severe and globally widespread infectious diseases affecting swine of all ages. PED has resulted in significant economic losses to the pig industry over the past three decades. It was first observed in Europe in 1971 [1]. Outbreaks occurred in Asian countries in 1982 and since then, PED has had a growing economic impact on the Asian swine industry [2]. In 2013, PEDV was first reported in the US and spread rapidly nationwide [3]. PED is characterized by severe diarrhea, vomiting, dehydration, and a mortality rate of up to 90% in suckling piglets [4].

Domestic pig farms in China have been using the epidemic diarrhea inactivated vaccine (KPEDV-9) since 2011. A more common method of immune prevention is to use infected sow feces and the intestines of infected piglets during pregnancy. The collected feces and intestines are artificially mixed and fed back (i.e., feedback method). However, the poor external biosecurity of pig farms makes the effectiveness of these preventive methods questionable because PED often occurs on farms that use vaccination or feedback methods. Therefore, research is now focused on alternative genetic-based methods, which first involves identifying specific host genes responsible for resistance to PED.

PEDV mainly infects and replicates in the villous enterocytes of the small intestine (duodenum, jejunum, and ileum) [5,6,7], which leads to villous atrophy and vacuolation, as well as a marked reduction in enzymatic activity. These changes interrupt the digestion and absorption of nutrients and electrolytes, leading to malabsorptive watery diarrhea in piglets [8]. The mechanism that leads to more serious disease and death following PEDV infection remains to be clearly defined [9]. In particular, receptor recognition of PEDV is still a subject of controversy. Over the last 10 years, similar to other coronaviruses, a number of studies have suggested that PEDV uses porcine aminopeptidase N (pAPN) as a cellular receptor [10]. Experiments showed biochemical interactions between PEDV S1 and APN expressed on the cell surface [11]. It can increase the susceptibility of these cells to PEDV infection when pAPN is expressed in non- tropic cell line (i.e, swine testis cells, ST cells) [12]. However, knockout of APN expression in PEDV-susceptible porcine and human cell lines (Huh7 and HeLa) confirmed that APN is not required for PEDV infection [13].

A few immune-related pathways, including the Janus kinase (JAK)/signal transducer and activator of transcription (STAT) signaling pathway, the NF-kappa B signaling pathway, and the PI3K/Akt/mTOR signaling pathway have been reported to play a role in responding to PEDV infection [14]. Importantly, the JAK-STAT signaling pathway regulates the adaptive and innate mechanisms related to mucosal immunity [15], and a previous study suggested that the JAK-STAT signaling pathway is activated after PEDV infection of host cells [16]. As an important component of the JAK-STAT signaling pathway, STAT proteins are reported to be activated in response to immunomodulation [17,18,19]. After activation, STAT proteins dimerize, translocate to the nucleus, and bind to the promoters of specific target genes, resulting in the regulation of expression of the target genes. To date, seven mammalian members of the STAT family have been identified, and have been suggested to play different physiological and biological roles.

Viral infection triggers changes in host gene expression patterns. In this study, to investigate the host response patterns to PEDV infection and to identify the key regulators involved in PEDV pathogenesis, the RNA-seq technique was applied to monitor global transcriptomic changes in the intestinal porcine epithelial cell line J2 (IPEC-J2). We report how signaling pathways respond to PEDV infection and suggest that STAT proteins act as key regulators when the host is under PEDV attack. Our findings offer insight into important host-dependent factors responsible for PEDV pathogenesis in vitro.

Results

Response of IPEC-J2 cells to PEDV strain CV777 infection

To detect the effects of PEDV strain CV777 on the growth phenotype of cells, healthy IPEC-J2 cells were infected with CV777 at a multiplicity of infection (MOI) of 1.0. As shown in Fig. 1a, the apoptosis rate increased with the increase in infection time. To further evaluate the response of IPEC-J2 cells to CV777 infection, we then detected the expression of the MX dynamin-like GTPase 1 (MX1) gene, a classic antiviral gene [20]. The results showed that MX1 mRNA expression was significantly increased after 48 h of CV777 infection (Fig. 1b). Next, the relative transcription level of the CV777 membrane protein (M) gene was also calculated. As shown in Fig. 1c, a significant increase in the expression level of M gene was detected 48 h after inoculation. These results suggested that it takes 48 h for CV777 to replicate extensively in IPEC-J2 cells and for IPEC-J2 cells to produce strong antiviral responses.

Fig. 1
figure 1

PEDV infection in IPEC-J2 cells. a Fluorescence-based cytotoxicity assay (using the LIVE/DEAD Viability/Cytotoxicity kit) of IPEC-J2 cells treated with PEDV strain CV777 at different time points. The green color indicates live cells stained by Calcein-AM and the red color indicates dead cells stained by propidium iodide. The green color was predominant in all viability assays, with only a few red (dead) cells appearing randomly. A significantly higher percentage of apoptotic cells (live (green)/dead (red)) are observed in the infected groups than in the controls (magnification, ×40; bars, 200 μm). b The mRNA expression level of MX1 gene in PEDV-infected IPEC-J2 cells as measured by qRT-PCR at different time points. c The mRNA expression level of M gene in PEDV-infected IPEC-J2 cells as measured by qRT-PCR at different time points. Data derive from three independent experiments and were analyzed by one-way ANOVA. Data are presented as the mean ± SD (n = 6). *** represents p < 0.001

RNA-seq analysis

To determine how global gene expression is differentially influenced by PEDV infection, RNA-seq was performed to compare gene expression between IPEC-J2 cells infected with PEDV for 48 h and mock-infected control cells. To avoid sample bias, six replicates for each group were collected and used for the construction of the RNA-seq libraries. The 12 libraries were sequenced on an Illumina HiSeq 2500 sequencer, and a total of 23.0–29.9 million paired raw reads were generated for each library. After removing low-quality reads and reads with adaptor sequences, we obtained approximately 22.9–29.9 million paired clean reads. Subsequently, the clean reads were aligned to the pig genome (Sus scrofa 11.1) and all samples had mapping ratios within the range of 91.86–92.99% (Table 1). Principal component analysis showed distinct separation between groups and confirmed the reproducibility of our biological samples (Fig. 2a). A sample correlation matrix based on gene expression levels further emphasized the high reproducibility and reliability of our experimental samples (Fig. 2b).

Table 1 Descriptive summary of data generated by RNA-Seq
Fig. 2
figure 2

Analysis of sample correlation. a PCA analysis based on the gene expression profile of each sample. b Heatmap shows the Pearson’s correlation of gene expression levels between samples

Identification of the host gene response to PEDV infection

A total of 18,573 pig genes was found to be expressed both in mock-infected cells (n = 6) and PEDV-infected cells (n = 6). On comparison, a total of 716 significantly upregulated and 138 significantly downregulated genes were identified in PEDV-infected cells compared with control cells (|log2 fold change| ≥ 1 and FDR ≤ 0.05; Fig. 3a, b; Supplementary Table 1), including many known antiviral genes, such as interferon-induced protein with tetratricopeptide repeats 1 (IFIT1) [21], MX1, myxovirus resistance 2 (MX2) [22], and tripartite motif containing 25 (TRIM25) [23]. The top five genes with the highest changes in the mRNA level were IFIT1, DExD/H-box helicase 58 (DDX58), radical S-adenosyl methionine domain containing 2 (RSAD2), MX2, and forkhead box S1 (FOXS1). To validate the quality of the RNA-seq data, five differentially expressed genes were random selected and their expression patterns were detected by qRT-PCR. The results showed that the expression patterns of the genes were consistent with the RNA-seq results (Fig. 3c), although the observed fold changes differed between the qRT-PCR and RNA-seq data, which may reflect differences in the sensitivity and specificity between qRT-PCR and high-throughput sequencing technology. Our findings suggested that the RNA-seq results were generally reliable.

Fig. 3
figure 3

Identification of differentially expressed genes. a Volcano plots of differentially expressed genes between the mock-infected and PEDV-infected groups. The red and green dots represent upregulated and downregulated genes, respectively, in the PEDV-infected group compared with the mock-infected group. b Heatmap showing the expression levels of the differentially expressed genes. Columns represent individual samples and rows indicate genes with significant expression differences between the two groups. c qRT-PCR validation of the differentially expressed genes. The left axis represents gene expression levels as determined by qRT-PCR, and the right axis represents the expression levels determined by RNA-seq in FPKM units. mRNA expression levels were normalized to the mRNA levels of the pig ACTB gene. Graphed data represent the mean ± SD, n = 6. * represents p < 0.05, ** represents p < 0.01, *** represents p < 0.001

Functional analysis of the differentially expressed genes

To better understand the functions of the differentially expressed genes, GO and KEGG pathway enrichment analyses were performed. The GO enrichment analysis showed that these genes were significantly enriched in 67 GO terms (p < 0.01, Supplementary Table 2). Among these, 19 (28.4%) GO terms were related to the host immune response and inflammatory response (Table 2), including the innate immune response (GO: 0045087), positive regulation of T cell proliferation (GO: 0042102), and regulation of the adaptive immune response (GO: 0002819). KEGG pathway enrichment analysis indicated that the differentially expressed genes were significantly enriched in 41 KEGG pathways (p < 0.01, Supplementary Table 3). Among these pathways, the NF-kappa B signaling pathway, Toll-like receptor (TLR) signaling pathway, JAK-STAT signaling pathway, and the intestinal immune network for IgA production signaling pathway were associated with the innate immune response and inflammatory response (Table 3). Interestingly, almost all differentially expressed genes that clustered in these immune-related GO and KEGG pathways were found to be upregulated after CV777 inoculation. For example, interferon regulatory factor 7 (IRF7), C-C motif chemokine ligand 5 (CCL5), STAT1, and interleukin 6 (IL6), which are clustered in the innate immune response (GO: 0045087), were all upregulated after PEDV infection. In the NF-kappa B signaling pathway, all differentially expressed genes, including CD40 molecule, TRIM25, and nuclear factor kappa B subunit 1 A (NFKB1A), were found to be upregulated. These results showed that the innate immune response and inflammatory response were activated during the course of CV777 infection.

Table 2 The GO terms were related to host immune response and inflammatory response
Table 3 The KEGG pathway were related to host immune response and inflammatory response

Transcription factor prediction among differentially expressed genes

The search for significantly overrepresented transcription factor binding sites in the promoter regions of the differentially expressed genes could be a powerful approach for finding key regulators of complex biological processes. Therefore, the putative promoter regions (2000 bp upstream of the transcription start site) of the differentially expressed genes were examined for the presence of transcription factor binding sites using the MEME tool. This analysis revealed several significantly enriched motifs (Supplementary Fig. 1), which were then annotated as the motif of transcription factors including zinc finger protein 384 (ZNF384), SRY-box transcription factor 10 (SOX10), STAT1 (a member of the STAT family), and recombination signal binding protein for immunoglobulin kappa J region (RBPJ), among others (Fig. 4a). In particular, five STAT members—STAT1, STAT2, STAT3, STAT4, and STAT5a—were identified as differentially expressed genes by RNA-seq (Fig. 4b). Subsequently, all of the STAT members and their known target genes, GBP1 [24, 25] and IFIT1 [26], were verified by qRT-PCR and the results were consistent with the RNA-seq data (Fig. 4b). To corroborate these findings, the Find Individual Motif Occurrences (FIMO) tool was used to identify putative STAT1 binding sites in these sequences using the count matrix motif of STAT obtained from the JASPAR database (ID: MA0137.3). As a control, another set of promoter sequences extracted from 854 randomly-selected genes were subjected to the same FIMO analysis. As a result, 504 (59.02%) promoter sequences for the differentially expressed genes were identified to have at least one STAT1 binding site. By comparison, only 112 (13.11%) of sequences among the control promoter sequences were identified as having a STAT1 binding site. The highly significant overrepresentation (p < 0.001, Pearson’s Chi-square test) of STAT1 transcription factors among the putative promoters of the differentially expressed genes suggested an important role for STAT factors during the PEDV infection process.

Fig. 4
figure 4

Prediction of potential transcription factors. a Enrichment of putative transcription factor binding motifs in the promoter regions of differentially expressed genes detected by MEME. The right panel shows the annotation of the enriched motifs as determined by MAST (ID, number, p value). b qRT-PCR validation of STAT expression levels and those of their known target genes. The left axis represents the gene expression levels determined by qRT-PCR, and the right axis represents the expression levels determined by RNA-seq in FPKM units. mRNA expression levels were normalized to the mRNA levels of the pig ACTB gene. Graphed data represent the mean ± SD, n = 6. * represents p < 0.05, ** represents p < 0.01, *** represents p < 0.001

Construction of a gene regulatory network between the STAT factors and the differentially expressed genes

On the basis of the assumption that the STAT protein family-mediated regulatory and signaling networks are representative of the infected interactome, STAT factors were used as seeds to construct a gene regulatory network. Figure 5a shows the nodes and interactions at the intersection of the network. GO and KEGG analyses showed that these genes were highly enriched in immune response-related functions (Fig. 5b). Additionally, topological analysis indicated that IL6, tumor necrosis factor (TNF), NFKB1A, MX1, and TLR2 are the principal hub genes, and all these genes are recognized as important immune response genes. Taken together, our results further implicate STAT members as key regulators when the host is under PEDV attack.

Fig. 5
figure 5

STAT transcription factors as key regulators during PEDV infection. a PPI network of the differentially expressed genes based on STRING analysis. The hub genes of the network are shown in yellow. b Top five pathways enriched in immune response-related functions. GO and KEGG analyses were conducted using a DAVID functional annotation program

Discussion

PEDV is a pathogen of interest to researchers because of its significant impact on the swine industry worldwide. Many previous studies have focused on viral isolation and molecular epidemiology surveys [27, 28]; however, it is also necessary to understand the molecular mechanisms involved in the host response to PEDV infection. Our study, which involved transcriptome analyses, revealed 854 significantly differentially expressed genes in the host. These differentially expressed genes were mainly enriched in the influenza A, TNF signaling pathway, inflammatory response, and other immune-related pathways. In particular, five members of the STAT family were significantly upregulated in PEDV-infected cells compared with mock-infected cells. STAT factors are considered to be downstream mediators of cytokine signaling. STAT proteins are capable of rapid and direct transduction of signals from the cell membrane to the nucleus, and they combine the functions of transduction agents and inducers of transcription. In the latent state, STAT proteins are inactive as monomers or unphosphorylated dimers, which are localized in the cytoplasm of unstimulated target cells [29,30,31,32]. As a non-negligible component of the JAK-STAT signaling pathway, STAT factors play an indispensable role in innate immunity to viral infection. Each member of the STAT family can be activated by multiple cytokines and their associated JAKs [33]. STAT1 plays a role in many important cytokine induction pathways and can upregulate many pro-inflammatory cytokines, whereas STAT2 is a co-factor in the type I IFN signaling pathway. Previous studies have shown that PEDV infection inhibits type I IFN induction [34]. Guo et al. demonstrated that PEDV infection promotes the degradation and interrupts the activation of STAT1 without inhibiting STAT1 transcription, while STAT2 remains uncleaved [35]. As the most important activator of STAT3, IL6 can directly act on target cells to affect STAT3 expression [36]. Yang et al. demonstrated that the direct interaction between PEDV S protein and epidermal growth factor receptor (EGFR) induces EGFR activation, and thus increases PEDV infection. They further demonstrated that activation of EGFR contributes to the enhancement and promotion of PEDV replication via JAK-STAT3 signaling pathways [37]. Our transcriptome data are consistent with previous studies. However, no significant change was observed in STAT5b and STAT6 expression after PEDV infection. Curiously, STAT5a and STAT5b, a pair of homologous genes that exhibit high similarity (~ 94%) in their coding sequences, showed different expression patterns after PEDV infection. STAT5a and STAT5b have multiple functions, including cell proliferation and differentiation [38], immunoregulation [39], a drug response, and metastasis [29]. STAT5a and STAT5b often perform similar functions in mediating regulatory signaling. However, recent studies have shown that STAT5a and STAT5b can also have distinct roles in regulating gene expression [40]. Lamba et al. reported the distinct and potentially opposing roles of STAT5a and STAT5b in the regulation of hepatic drug response genes [41]. STAT5a is expressed at a much lower level in the liver than STAT5b [42]. Jennifer et al. reported the differentially regulated expression of FOXP3 and IL-2R in STAT5b knockdown human primary T cells and the downregulated expression of Bcl-X only in STAT5a knockdown human primary T cells [43]. In the present study, our findings revealed for the first time that STAT5a, but not STAT5b, regulated the expression of key factors involved in the immune response after virus infection.

Previous studies provided transcriptomic and proteomic data for PEDV-infected Vero cells that offered the first valuable insight to better understand the host response to PEDV infection [16, 44, 45]. However, these reports contained some contradictory results. For example, nearly all of the detected IRFs were reported to be downregulated after PEDV infection in Vero E6 cells [44], whereas our results showed that IRF1, FIRF2, IRF7, and IRF9 were significantly upregulated after PEDV infection in IPEC-J2 cells, and other detected IRF members were not significantly changed. These conflicting results are most likely due to the different cell lines or the different detection times after infection used between studies. With this in mind, we believe that more integrated and detailed studies of time-resolved or even single-cell-resolved gene expression changes will aid in understanding the host gene response to PEDV infection. A recent study performed global mapping of H3K4 trimethylation and transcriptomic analyses in the PEDV-infected jejunum of piglets compared with healthy piglets, using chromatin immunoprecipitation sequencing and RNA-seq, and provided novel insights into our understanding of the host response to PEDV infection [14]. RSAD2, an important antiviral gene [46], has attracted extensive attention recently. Jiang et al. confirmed that RSAD2 is necessary for mouse dendritic cell (mDC) maturation via the IRF7-mediated signaling pathway, and mDCs were shown to lose their antitumor efficacy in a RSAD2-knockdown lung metastasis mouse model [47]. IRFs are involved in antiviral defense and immune regulation [48, 49]. In particular, IRF3 and IRF7 are direct transducers of virus-mediated signaling in the induction of type I IFN [50]. Our results showed that the mRNA expression of RASD2 and IRF7 was significantly upregulated, while the mRNA expression of IRF3 did not change after PEDV infection. This result suggested that IRF7, but not IRF3, might be regulated by RSAD2 to initiate the antiviral response.

Viruses have the ability to manipulate host cell signal transduction pathways to promote their own survival [51]. In our study, significantly enriched pathways in response to PEDV infection included the protein kinase B signaling (Akt) pathway, which modulates the cell cycle, apoptosis, and differentiation. It is well established that the regulation of Akt activity is mainly dependent upon PI3K activity. Unexpectedly, phosphoinositide-3-kinase adaptor protein 1 (PIK3AP1), a member of the vital PIKs, was significantly upregulated during PEDV infection at the transcriptional level. The mechanism by which PEDV manipulates the PIK/Akt pathway to achieve maximum virus replication is worth studying.

Conclusion

In this study, changes in a series of immune and inflammatory response-related genes in PEDV CV777-infected IPEC-J2 were observed by RNA-seq analysis. Among the total gene set, 59.02% of fragments were identified as possessing at least one putative binding site for STAT transcription factors. Moreover, a STAT member-induced regulatory network, including 53 genes that were differentially expressed during the process of PEDV infection, was constructed. These results will be useful for exploring and further understanding host responses to PEDV infection in pigs.

Methods

Cells and virus

IPEC-J2 cells were grown in antibiotic-free Dulbecco’s modified Eagle’s medium (DMEM)/F12 (Gibco, Rockville, MD, USA) supplemented with 5% fetal bovine serum (Gibco), insulin (5 μg/mL, Gibco), transferrin (5 μg/mL, Gibco), selenium (5 ng/mL, Gibco), and epidermal growth factor (5 ng/mL, Sigma-Aldrich, St. Louis., MO, USA). Cells were maintained in 100 mL flasks at 37 °C in an atmosphere of 5% CO2. The medium was replaced every 2 days.

CV777 strain of porcine epidemic diarrhea virus was provided by the College of Animal Medicine of China Agricultural University. The IPEC-J2 line was initially obtained from the American Type Culture Collection (ATCC, Manassas, VA, USA).

PEDV infection of IPEC-J2 cells

Culture medium was removed when the monolayer of IPEC-J2 cells had grown to 80% confluence in 10 cm2 cell culture dishes. After washing gently with PBS, cells were incubated in serum-free DMEM and inoculated with CV777 at a multiplicity of infection of 1. The serum-free DMEM was replaced with culture medium after 1 h of adsorption. Cells were harvested after another 12, 24, 36, and 48 h of culture, and were infected with PEDV. Cells incubated in serum-free DMEM without inoculation of CV777 served as the mock-infected control group.

Cell apoptosis detection

Cell apoptosis was detected using the TransDetect Cell LIVE/DEAD Viability/Cytotoxicity Detection Kit (Trans, Beijing, China) according to the manufacturer’s protocol. Briefly, culture medium was removed, and cells were washed once with PBS. Cells were then incubated in a mixture of 1 μL of 1000× Calcein-AM solution, 1 μL of 1000× propidium iodide (PI) working solution, and 500 μL of PBS for 15 min at 37 °C in the dark. Cells were then immediately observed using a fluorescence microscope (Nikon, Tokyo, Japan). Live cells were marked with green fluorescence and apoptotic cells were marked with red fluorescence.

RNA-seq

RNA-seq was performed to compare gene expression profiles between PEDV-infected IPEC-J2 cells and uninfected IPEC-J2 cells over a 48 h time period. Six biological repeats for each group were performed to prepare the RNA-seq library. The RNeasy mini kit (Qiagen, Hamburg, Germany) was used to isolate total RNA from the 12 samples according to the manufacturer’s instructions. Total RNA was purified using the RNeasy micro kit (Qiagen) and RNase-Free DNase Set (Qiagen). The RNA concentration was determined using the NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific Inc., USA). RNA integrity was checked using the RNA Integrity Number (RIN) generated by an Agilent Bioanalyzer 2100. RNA samples with RIN ≥ 7 were used for library construction. Sequencing libraries were constructed using the TruSeq®RNA Sample Preparation Kit v2 (Illumina, San Diego, CA, USA) following the manufacturer’s protocol. Briefly, RNA was initially purified using poly-T oligo-attached magnetic beads. The purified RNA was fragmented into 100–400 bp small pieces using divalent cations at 94 °C for 8 min. The double-stranded cDNA was synthesized by priming with random hexamers. A single nucleotide A (adenine) was added (A-tailing) to the 3′ end of the end-repaired DNA fragments. The repaired fragments were ligated to sequencing adapters. The products were purified by agarose gel electrophoresis and were used as template in PCR. The minimal number of PCR cycles (15 cycles) was used to avoid skewing the representation of the library. The final cDNA libraries were quantified using a Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and validated using an Agilent Bioanalyzer 2100 (Agilent Technologies). Finally, 12 libraries were sequenced on the Illumina HiSeq 2500 sequencer (Illumina). Transcriptome sequencing was entrusted to Beijing Boao Jingdian Biological Co., Ltd.

Sequence bioinformatics

Raw data were cleaned using BBMap (https://sourceforge.net/projects/bbmap/) to remove sequencing linkers, primer sequences, and low-quality reads. Then, the quality-controlled sequencing data were aligned to the downloaded pig reference genome using TopHat2 software [52]. The pig reference genome (Sus scrofa 11.1) and the corresponding annotation file were downloaded from the Ensembl database. Then, Cufflinks software was employed to quantify the fragment per kilobase of exon model per million mapped reads (FPKM) values for each gene. The Cuffdiff software was used to screen the differentially expressed genes [53] . The significance threshold for the differential expression was FDR < 0.05 and a |log2 fold change | > 1.

GO and KEGG pathway enrichment

Functional enrichment of the gene module was analyzed using the web-based tools in DAVID (v6.8; https://david.ncifcrf.gov/) to identify enriched GO terms and KEGG pathways. Ensembl gene IDs were submitted to the Gene Functional Classification Tool, then biological process, cellular component, molecular function, and KEGG pathway were selected to perform the enrichment. GO terms and KEGG pathways with p < 0.05 were defined as significantly enriched.

Transcription factor binding motif analysis

We interrogated enrichments of transcription factors in the promoters of differentially expressed genes using the MEME Suite (http://meme-suite.org/index.html). First, the putative promoter regions (2000 bp upstream of the transcription start site) of the differentially expressed genes were extracted using the BioMart program in the Ensembl database. Then, these promoter sequences were employed to discover enriched motifs using the MEME tool under the default parameters. The discover program was run repeatedly until the five most enriched motifs were discovered. Subsequently, the discovered motifs were compared against the JASPAR CORE database using the Tomtom tool, resulting in exact motif annotations. Additionally, the FIMO tool was employed to scan the set of promoter sequences for individual matches to specific motifs obtained from the JASPAR CORE database. FIMO analysis was performed with stringent criteria including a p value <5E˗4 and a maximum of two mismatched residues.

Construction of a gene regulatory network

The STRING (Search Tool for the Retrieval of Interacting Genes, http://string-db.org/, version 11.0) database, which provides both predicted and experimental interaction information for proteins [54], was used to construct the gene regulatory network. Interactions with a combined score > 0.4 were considered statistically significant. We then exported the results to Cytoscape v3.7.2 to visualize the gene regulatory network. In the network, each node represents a protein and each edge represents an interaction between two proteins.

Quantitative real-time PCR

One microgram of RNA from each sample was used to synthesize cDNA using the PrimeScriptâ„¢ RT reagent kit with gDNA Eraser (TaKaRa Biotech, Dalian, China) according to the manufacturer’s instructions. qRT-PCR analysis was carried out using SYBR Green master mix (Roche, USA). Each reaction contained 5 μL of SYBR Green PCR Master Mix (TaKaRa), 0.5 μL of forward primer, 0.5 μL of reverse primer, 0.5 μL of cDNA, and 3.5 μL of RNase-free water. The qRT-PCR program was as follows: 95 °C for 5 min; followed by 40 PCR cycles at 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 20 s; and then a further 10-min extension at 72 °C. Each individual sample for qRT-PCR was run in triplicate. The gene expression levels were measured using the 2˗△△Ct method. The pig ACTB gene, a housekeeping gene, was used as an internal control. All primers were designed using the NCBI Primers-BLAST online programs (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) and their sequences are listed in Supplementary Table 4.

Statistical analysis

Statistical analyses of the qRT-PCR results were carried out using SPSS 20.0 (IBM Corp., Armonk, NY, USA). The gene expression levels are presented as the mean ± SD and differences between groups were evaluated using one-way ANOVA (Dunnett’s t-test) and a two-tailed Student’s t-test.

Availability of data and materials

The datasets generated during the current study are available in the NCBI Sequence Read Archive under the accession number PRJNA599091.

Abbreviations

PED:

Porcine epidemic diarrhea

PEDV:

Porcine epidemic diarrhea virus

IPEC-J2:

Intestinal porcine epithelial cell lines J2

STAT:

Signal transducer and activator of transcription

JAK:

Janus kinase

pAPN:

Porcine aminopeptidase N

MX1:

MX dynamin like GTPase 1

IFIT1:

Interferon-induced protein with tetratricopeptide repeats 1

MX2:

Myxovirus resistance 2

TRIM25:

Tripartite motif containing 25

DDX58:

DExD/H-box helicase 58

RSAD2:

Radical S-adenosyl methionine domain containing 2

FOXS1:

Forkhead box S1

IRF7:

Interferon regulatory factor 7

CCL5:

C-C motif chemokine ligand 5

IL6:

Interleukin 6

NFKB1A:

Nuclear factor kappa B subunit 1 A

ZNF384:

Zinc finger protein 384

SOX10:

SRY-box transcription factor 10

RBPJ:

Recombination signal binding protein for immunoglobulin kappa J region

FIMO:

Find Individual Motif Occurences

EGFR:

Epidermal growth factor receptor

References

  1. Wood EN. An apparently new syndrome of porcine epidemic diarrhoea. Veterinary Record. 1977;100(12):243–4.

    Article  CAS  Google Scholar 

  2. Song D, Park B. Porcine epidemic diarrhoea virus: a comprehensive review of molecular epidemiology, diagnosis, and vaccines. Virus Genes. 2012;44(2):167–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Cima G. Fighting a deadly pig disease. Industry, veterinarians trying to contain PED virus, new to the US. J Am Vet Med Assoc. 2013;243(4):469–70.

    PubMed  Google Scholar 

  4. Pensaert MB, Bouck PD. A new coronavirus-like particle associated with diarrhea in swine. Arch Virol. 1978;58(3):243–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Debouck P. The pathogenesis of an enteric infection in pigs. Experimentally induced by the coronavirus-like agent, CV 777. Vet Microbiol. 1981;6(2):157–65.

    Article  PubMed Central  Google Scholar 

  6. Ducatelle R, Coussement W, Charlier G, Debouck P, Hoorens J. Three-dimensional sequential study of the intestinal surface in experimental porcine CV 777 coronavirus enteritis. Zentralblatt Für Veterinärmedizinreihe Bjournal of Veterinary Medicine. 1981;28(6):483–93.

    Article  CAS  Google Scholar 

  7. Ducatelle R, Coussement W, Pensaert MB, Debouck P, Hoorens J. In vivo morphogenesis of a new porcine enteric coronavirus, CV 777. Arch Virol. 1981;68(1):35–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Ducatelle R, Coussement W, Debouck P, Hoorens J. Pathology of experimental CV777 coronavirus enteritis in piglets. II. Electron microscopic study. Vet Pathol. 1982;19(1):57–66.

    Article  CAS  PubMed  Google Scholar 

  9. Sun D, Wang X, Wei S, Chen J, Feng L. Epidemiology and vaccine of porcine epidemic diarrhea virus in China: a mini-review. J Veterinary Med Sci. 2016;78(3):355–63.

  10. Li BX, Ge JW, Li YJ. Porcine aminopeptidase N is a functional receptor for the PEDV coronavirus. Virology. 2007;365(1):166–72.

    Article  CAS  PubMed  Google Scholar 

  11. Deng F, Ye G, Liu Q, Navid MT, Zhong X, Li Y, et al. Identification and comparison of receptor binding characteristics of the spike protein of two porcine epidemic diarrhea virus strains. Viruses. 2016;8(3):55.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Nam E, Lee C. Contribution of the porcine aminopeptidase N (CD13) receptor density to porcine epidemic diarrhea virus infection. Vet Microbiol. 2010;144(1–2):41–50.

    Article  CAS  PubMed  Google Scholar 

  13. Li W, Luo R, He Q, van Kuppeveld FJM, Rottier PJM, Bosch BJ. Aminopeptidase N is not required for porcine epidemic diarrhea virus cell entry. Virus Res. 2017;235:6–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Wang H, Yang L, Qu H, Feng H, Wu S, Bao W. Global mapping of H3K4 Trimethylation (H3K4me3) and Transcriptome analysis reveal genes involved in the response to epidemic diarrhea virus infections in pigs. Animals. 2019;9(8):523.

    Article  CAS  PubMed Central  Google Scholar 

  15. Heneghan AF, Pierre JF, Kudsk KA. JAK-STAT and intestinal mucosal immunology. Jak-stat. 2013;2(4):e25530.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Lin H, Li B, Chen L, Ma Z, He K, Fan H. Differential protein analysis of IPEC-J2 cells infected with porcine epidemic diarrhea virus pandemic and classical strains elucidates the pathogenesis of infection. J Proteome Res. 2017;16(6):2113–20.

    Article  PubMed  CAS  Google Scholar 

  17. Sadowski HB, Shuai K, Darnell JE Jr, Gilman MZ. A common nuclear signal transduction pathway activated by growth factor and cytokine receptors. Science. 1993;261(5129):1739–44.

    Article  CAS  PubMed  Google Scholar 

  18. Wegenka UM, Lütticken C, Buschmann J, Yuan J, Lottspeich F, Müller-Esterl W, et al. The interleukin-6-activated acute-phase response factor is antigenically and functionally related to members of the signal transducer and activator of transcription (STAT) family. Mol Cell Biol. 1994;14(5):3186–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Darnell JE Jr, Kerr IM, Stark GR. Jak-STAT pathways and transcriptional activation in response to IFNs and other extracellular signaling proteins. Science. 1994;264(5164):1415–21.

    Article  CAS  PubMed  Google Scholar 

  20. Pillai PS, Molony RD, Martinod K, Dong H, Pang IK, Tal MC, et al. Mx1 reveals innate pathways to antiviral resistance and lethal influenza disease. Science. 2016;352(6284):463–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Barik S. In silico structure analysis of alphaviral RNA genomes shows diversity in the evasion of IFIT1-mediated innate immunity. J Biosci. 2019;44(4):79.

    Article  PubMed  CAS  Google Scholar 

  22. Betancor G, Dicks MDJ, Jimenez-Guardeno JM, Ali NH, Apolonia L, Malim MH. The GTPase Domain of MX2 Interacts with the HIV-1 Capsid, Enabling Its Short Isoform to Moderate Antiviral Restriction. Cell Rep. 2019;29(7):1923–33 e3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Chen ST, Chen L, Lin DS, Chen SY, Tsao YP, Guo H, et al. NLRP12 Regulates Anti-viral RIG-I Activation via Interaction with TRIM25. Cell Host Microbe. 2019;25(4):602–16 e7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Müller P, Pugazhendhi D, Zeidler MP. Modulation of human JAK-STAT pathway signaling by functionally conserved regulators. Jak-stat. 2012;1(1):34–43.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Ehret GB, Reichenbach P, Schindler U, Horvath CM, Fritz S, Nabholz M, et al. DNA binding specificity of different STAT proteins. Comparison of in vitro specificity with natural target sites. J Biol Chem. 2001;276(9):6675–88.

    Article  CAS  PubMed  Google Scholar 

  26. Zhou X, Michal JJ, Zhang L, Ding B, Lunney JK, Liu B, et al. Interferon induced IFIT family genes in host antiviral defense. Int J Biol Sci. 2013;9(2):200–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Pan Y, Tian X, Li W, Zhou Q, Wang D, Bi Y, et al. Isolation and characterization of a variant porcine epidemic diarrhea virus in China. Virol J. 2012;9:195.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Chen Q, Li G, Stasko J, Thomas JT, Stensland WR, Pillatzki AE, et al. Isolation and characterization of porcine epidemic diarrhea viruses associated with the 2013 disease outbreak among swine in the United States. J Clin Microbiol. 2014;52(1):234–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Xiong H, Su W-Y, Liang Q-C, Zhang Z-G, Chen H-M, Du W, et al. Inhibition of STAT5 induces G1 cell cycle arrest and reduces tumor cell invasion in human colorectal cancer cells. Lab Investig. 2009;89(6):717–25.

    Article  CAS  PubMed  Google Scholar 

  30. Lim CP, Cao X. Structure, function, and regulation of STAT proteins. Mol BioSyst. 2006;2(11):536–50.

    Article  CAS  PubMed  Google Scholar 

  31. Haan S, Kortylewski M, Behrmann I, Muller-Esterl W, Heinrich PC, Schaper F. Cytoplasmic STAT proteins associate prior to activation. Biochemical J. 2000;345(3):417–21.

    Article  CAS  Google Scholar 

  32. Braunstein J, Brutsaert S, Olson R, Schindler C. STATs Dimerize in the absence of phosphorylation. J Biol Chem. 2003;278(36):34133–40.

    Article  CAS  PubMed  Google Scholar 

  33. O'Shea JJ, Schwartz DM, Villarino AV, Gadina M, McInnes IB, Laurence A. The JAK-STAT pathway: impact on human disease and therapeutic intervention. Annu Rev Med. 2015;66:311–28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Zhang Q, Shi K, Yoo D. Suppression of type I interferon production by porcine epidemic diarrhea virus and degradation of CREB-binding protein by nsp1. Virology. 2016;489:252–68.

    Article  CAS  PubMed  Google Scholar 

  35. Guo L, Luo X, Li R, Xu Y, Zhang J, Ge J, et al. Porcine epidemic diarrhea virus infection inhibits interferon signaling by targeted degradation of STAT1. J Virol. 2016;90(18):8281–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Johnson DE, RA OK, JR G. Targeting the IL-6/JAK/STAT3 signalling axis in cancer. Nat Rev Clin Oncol. 2018;15(4):234.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Yang L, Xu J, Guo L, Guo T, Zhang L, Feng L et al. Porcine Epidemic Diarrhea Virus-Induced Epidermal Growth Factor Receptor Activation Impairs the Antiviral Activity of Type I Interferon. J Virol. 2018;92(8):e02095–17.

  38. Jin HP, Darvin P, Lim EJ, Joung YH, Hong DY, Park EU, et al. Hwanggeumchal sorghum induces cell cycle arrest, and suppresses tumor growth and metastasis through Jak2/STAT pathways in breast Cancer Xenografts. PLoS One. 2012;7(7):e40531.

    Article  CAS  Google Scholar 

  39. Heltemes-Harris LM, Farrar MA. The role of STAT5 in lymphocyte development and transformation. Curr Opin Immunol. 2012;24(2):146–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Lamba V, Ghodke Y, Guan W, Tracy TS. microRNA-34a is associated with expression of key hepatic transcription factors and cytochromes P450. Biochem Biophys Res Commun. 2014;445(2):404–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Lamba V, Jia B, Liang F. STAT5A and STAT5B have opposite correlations with drug response gene expression. Biochemical & Biophysical Research Communications. 2016;479(2):117–24.

    Article  CAS  Google Scholar 

  42. Park SH, Liu X, Hennighausen L, Davey HW, Waxman DJ. Distinctive roles of STAT5a and STAT5b in sexual dimorphism of hepatic P450 gene expression. Impact of STAT5a gene disruption. J Biol Chem. 1999;274(11):7421–30.

    Article  CAS  PubMed  Google Scholar 

  43. Jenks JA, Seki S, Kanai T, Huang J, Morgan AA, Scalco RC, et al. Differentiating the roles of STAT5B and STAT5A in human CD4+ T cells. Clin Immunol. 2013;148(2):227–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Zhang H, Liu Q, Su W, Wang J, Sun Y, Zhang J, et al. Genome-wide analysis of differentially expressed genes and the modulation of PEDV infection in Vero E6 cells. Microb Pathog. 2018;117:247–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Zeng, H Z, Z D, R L, K A, L L et al. proteome analysis of porcine epidemic diarrhea virus (PEDV)-infected Vero cells. Proteomics. 2015;15(11):1819–1828.

  46. Honarmand EK. A unifying view of the broad-spectrum antiviral activity of RSAD2 (viperin) based on its radical-SAM chemistry. Metallomics. 2018;10(4):539–52.

    Article  Google Scholar 

  47. Ji-Su J, Jun-Ho L, Nam-Chul J, So-Yeon C, Soo-Yeoun P, Ji-Young Y, et al. Rsad2 is necessary for mouse dendritic cell maturation via the IRF7-mediated signaling pathway. Cell Death Dis. 2018;9(8):823.

    Article  CAS  Google Scholar 

  48. Honda K, Taniguchi T. IRFs: master regulators of signalling by toll-like receptors and cytosolic pattern-recognition receptors. Nat Rev Immunol. 2006;6(9):644–58.

    Article  CAS  PubMed  Google Scholar 

  49. Nguyen H, Hiscott J, Pitha PM. Cytokine Growth Factor Rev. 1997;8(4):293–312.

    Article  CAS  PubMed  Google Scholar 

  50. Marie I. Differential viral induction of distinct interferon-alpha genes by positive feedback through interferon regulatory factor-7. EMBO J. 1998;17(22):6660–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Bonjardim CA. Viral exploitation of the MEK/ERK pathway - a tale of vaccinia virus and other viruses. Virology. 2017;507:267–75.

    Article  CAS  PubMed  Google Scholar 

  52. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Ghosh S, Chan CK. Analysis of RNA-Seq Data Using TopHat and Cufflinks. 2016;1374:339–61.

    CAS  Google Scholar 

  54. von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31(1):258–61.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by the Shenzhen Key Technology Projects (JSGG20180507182028625).

Author information

Authors and Affiliations

Authors

Contributions

JL was responsible for study conception and design; ZH, YL and HD analyzed and interpreted the patient data. ZH drafted the manuscript. JL, JR, XZ and KW revised of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Jianfeng Liu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplementary Figure 1.

Motifs discovered using the MEME tool.

Additional file 2: Supplementary Table 1.

Information on each significantly differentially expressed gene (p < 0.05) identified in PEDV-infected cells.

Additional file 3: Supplementary Table 2.

GO terms significantly enriched among the differentially expressed genes in PEDV-infected cells.

Additional file 4: Supplementary Table 3.

KEGG pathways significantly enriched among the differentially expressed genes in PEDV-infected cells.

Additional file 5: Supplementary Table 4.

List of qRT-PCR primers.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hu, Z., Li, Y., Du, H. et al. Transcriptome analysis reveals modulation of the STAT family in PEDV-infected IPEC-J2 cells. BMC Genomics 21, 891 (2020). https://doi.org/10.1186/s12864-020-07306-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-020-07306-2

Keywords