RNA-Seq-Based Whole Transcriptome Analysis of IPEC-J2 Cells During Swine Acute Diarrhea Syndrome Coronavirus Infection

The new emergence of swine acute diarrhea syndrome coronavirus (SADS-CoV) has resulted in high mortality in suckling pigs in China. To date, the transcriptional expression of host cells during SADS-CoV infection has not been documented. In this study, by means of RNA-Seq technology, we investigated the whole genomic expression profiles of intestinal porcine epithelial cells (IPEC-J2) infected with a SADS-CoV strain SADS-CoV-CH-FJWT-2018. A total of 24,676 genes were identified: 23,677 were known genes, and 999 were novel genes. A total of 1,897 differentially expressed genes (DEGs) were identified between SADS-CoV-infected and uninfected cells at 6, 24, and 48 h post infection (hpi). Of these, 1,260 genes were upregulated and 637 downregulated. A Gene Ontology enrichment analysis revealed that DEGs in samples from 6, 24, and 48 hpi were enriched in 79, 383, and 233 GO terms, respectively, which were mainly involved in immune system process, response to stimulus, signal transduction, and cytokine–cytokine receptor interactions. The 1,897 DEGs were mapped to 109 KEGG Ontology (KO) pathways classified into four main categories. Most of the DEGs annotated in the KEGG pathways were related to the immune system, infectious viral disease, and signal transduction. The mRNA of porcine serum amyloid A-3 protein (SAA3), an acute phase response protein, was significantly upregulated during the infection. Over-expressed SAA3 in IPEC-J2 cells drastically inhibited the replication of SADS-CoV, while under-expressed SAA3 promoted virus replication. To our knowledge, this is the first report on the profiles of gene expression of IPEC-J2 cells infected by SADS-CoV by means of RNA-Seq technology. Our results indicate that SADS-CoV infection significantly modified the host cell gene expression patterns, and the host cells responded in highly specific manners, including immune response, signal and cytokine transduction, and antiviral response. The findings provide important insights into the transcriptome of IPEC-J2 in SADS-CoV infection.

Previous pathogenesis studies on coronaviruses have revealed that there are extensive and complex interactions between viruses and hosts, which is of import in terms of disease prevention and control caused by these pathogens. Host cell responses to virus infection involve complex interactions between cellular and viral networks (10)(11)(12). Most viruses attempt to change host cellular processes and organism systems to improve the efficiency of virus infection, and, on the other hand, the cells employ multiple mechanisms in responses to generating an antiviral state (13,14). As reported, CoV infection can cause alterations in the transcription and translation patterns, cell cycle, cytoskeleton, and apoptosis pathways of host cells (15)(16)(17)(18)(19)(20)(21)(22). With traditional approaches, it is difficult to explore the intricate and mass interactions between viruses and hosts/cells. Next-generation' sequencing (NGS) technology has provided a powerful tool for the studies of emerging and re-emerging human and animal pathogens/diseases. NGS technology, including whole genome sequencing, RNA-Seq, and single cell sequencing, has been extensively applied in genome sequencing, gene expression profiling analysis, and pathogen detection (15,(23)(24)(25).
Up to date, the study on mechanisms of infection and pathogenesis of SADS-CoV is limited. A study reported that SADS-CoV infection could antagonize interferon production and lead to immune evasion (26). However, the underneath mechanisms remain roughly unknown. There is no report on the transcriptome profile of intestinal epithelial cells during SADS-CoV infection. Thus, to address the global gene expressions profile of intestinal epithelial cells in SADS-CoV infection, a porcine intestinal epithelial cell line of IPEC-J2 was used as a model. The whole transcriptomics of the cells at 6, 24, and 48 h post infection (hpi) of SADS-CoV was determined via RNA-Seq technology. Subsequently, differentially expressed genes (DEGs) were screened, and the gene functions, signaling pathways associated with viral infection, and pathogenesis were analyzed. This study would increase our knowledge on the transcriptomics landscape of SADS-CoV infected small intestinal cells and shed light on future explorations of the mechanisms of SADS-CoV infection.

Growth Kinetics of SADS-CoV in IPEC-J2 Cells
To determine the growth kinetics of the isolated SADS-CoV strain SADS-COV-CH-FJWT-2018, confluent IPEC-J2 cells in six-well plates were washed twice with D-Hanks and then inoculated with 500 µl of viral supernatants containing SADS-COV-CH-FJWT-2018 (P10) at an MOI of 1. After 2 h of incubation at 37 • C with 5% CO 2 , the six-well plates were washed with D-Hanks to remove unabsorbed viruses and added 2 ml of DMEM containing 10 µg/ml trypsin. Afterwards, the six-well plates were incubated at 37 • C with 5% CO 2 , and cell culture supernatants were then collected at different time intervals (0, 6, 12, 18, and up to 72 hpi). The virus titer for each time point in each sample collected was determined by 50% tissue culture infective dose (TCID 50 ) (29) and an SYBR Green real-time RT-PCR assay previously established in our laboratory. Samples harvested at each time point were independently repeated three times, and the mean value, and standard deviation (SD) were calculated.

RNA Extraction, cDNA Library Construction, and RNA-Seq
For RNA-Seq, nine samples from IPEC-J2 cells infected with SADS-CoV at 6, 24, and 48 hpi with three samples from each time-point (group named as A_6h, A_24h, and A_48h, respectively), and nine mock-infected negative controls with corresponding time points (group named as B_6h, B_24h, and B_48h, respectively), were collected. Firstly, the cell culture supernatants were removed, and the cells were washed twice with 0.01 M sterile phosphate-buffered saline (pH 7.2). The cells were then lysed with 1 ml Trizol Reagent (TaKaRa, Japan) for RNA extraction according to the manufacturer's instructions. The concentrations of extracted RNA were determined by using a NanoDrop 2,000 spectrophotometer (Thermo Scientific, USA), and the integrity of purified RNA was evaluated via a RNA integrity number (RIN) by an Agilent Bioanalyzer 2,100 system (Agilent Technologies, USA). Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer's protocols. Briefly, mRNA was firstly purified from total RNA using poly-T oligo-attached magnetic beads, and then digested into short fragments by adding fragmentation reagents. Subsequently, first strand cDNA was synthesized using random N6 primer and M-MLV Reverse Transcriptase. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3 ′ ends of DNA fragments, adaptors with a hairpin loop structure were ligated to prepare for hybridization (NEBNest Adaptor Kit). Then, cDNA fragments were separated by agarose gel electrophoresis, and the fragments of 250-300 bp in length were selected and purified with AMPure XP system (Beckman Coulter, USA). Then, a 3 µl measurement of USER Enzyme (NEB, USA) was incubated with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR. Then PCR was performed to generate a cDNA library with Phusion High-Fidelity DNA polymerase, universal PCR primers, and Index (X) Primers (NEB lab, USA). Finally, PCR products were purified (AMPure XP system), and library quality was assessed on the Agilent Bioanalyzer 2,100 system. The libraries were sequenced on Illumina HiSeq 2,500 platform using the paired-end technology by Novogene Co., Ltd (Beijing, China).

Bioinformatic Analysis of RNA-Seq Data
The raw reads were qualified by removing reads containing adapter, poly-N and low-quality. The plus reads were clean reads, and all of the clean reads were then separated according to the barcodes. The adapter and barcode sequences were trimmed. The trimmed clean reads were assembled and mapped to the pig genome (WASHUC2.69) in Ensemble using Hisat2 (v2.0.5), and the fragments per kilobase million (FPKM) of each gene were then calculated based on the length of the gene and read counts mapped to genes. A negative binomial distribution-based model in DESeq2 was used to determine the differential expressed genes (DEGs) among the treatments. The fold change of genes between SADS-CoV infected and mock-infected samples were calculated under log 2 |FoldChange| (log 2 |FC|) ≥ 1 and a false discovery rate (FDR) adjusted p (p adj ) < 0.05 based on a method described previously (30). The numbers of novel and known genes were statistically calculated. To predict the major biological and molecular functions of these DGEs, Gene Ontology (GO) enrichment analysis (http://www.geneontology.org) was carried out and visualized by the cluster Profiler R package (31). All of the DEGs generated were functionally analyzed against KEGG (32). The significance of all GO and KEGG terms was corrected by controlling the p adj -value.

Verification of Differential Expressed Genes by Real-Time qRT-PCR
To validate the accuracy of the RNA-Seq results, 20 biologically related DEGs (IFIT1, IFIT2, IL-6, etc.) and a housekeeping gene beta actin were randomly selected for real-time qRT-PCR validation (primer information is listed in Additional file 1: Table S1). The relative expression values were normalized, with the beta actin gene serving as an internal control. After amplification, the relative fold change of the differentially expressed genes was calculated through the 2 − CT algorithm. and normal IPEC-J2 cells was measured by TCID 50 assays. SPSS software was used for statistical analysis (IBM, USA). Significant differences between groups were evaluated by using the Student's t-test. A threshold of p < 0.05 was regarded as a significant difference between two groups.

SADS-CoV Infection in IPEC-J2 Cells and Vero-81 Cells
Vero 81 and IPEC-J2 cell lines were routinely employed to propagate the cell-adapted SADS-COV-CH-FJWT-2018. In Vero-81 cells, the CPE induced by SADS-CoV was characterized by cell fusion, while in IPEC-J2 cells, the CPE displayed as enlarged, rounded, and condensed granular particles ( Figure 1A). To confirm the characteristics of SADS-CoV-CH-FJWT-2018 growth on IPEC-J2 cells, a growth curve was determined by quantifying the genomic RNA copies at different time intervals (0, 6, 12, 18, and up to 72 hpi) ( Figure 1B). From 0 to 24 h, the viral titers were increased and reached a highest titer of 10 8.706 copies/µl, and then decreased to 10 7.006 copies/µl from 36 to 72 h. The results demonstrated that IPEC-J2 cells are susceptible and permissive to SADS-CoV-CH-FJWT-2018, and the cell line is a good candidate for SADS-CoV research.
Expression Profiles of IPEC-J2 Cells Infected With SADS-COV-CH-FJWT-2018 The raw reads of RNA-Seq were quality controlled to ensure that all data were met the criteria for the whole transcriptomic analysis. A total of 132.80 Gb qualified bases were obtained with an average of 7.38 ± 0.70 Gb in each sample (Additional file 2: Table S2). The clean reads were then mapped to the reference pig genome (WASHUC2.69) in Ensemble using Hisat2 (v2.0.5), with an average of 94.84% qualified reads mapped reference sequences. A total of 24,676 genes were identified in which 23,677 were known genes and 999 novel genes (Additional file 3: Table S3).
Of these known genes, 17,745 were protein coding genes, 4,096 were long non-coding RNA, 186 were miRNA, and 1,237 were pseudogenes or transcript pseudogenes.

Differentially Expressed Genes Analyses
To identify the DEGs in response to SADS-CoV infection in IPEC-J2 cells, all the gene numbers were homogenized by an algorithm of Reads Per kb per Million reads (RPKM), and then DEGs were generated by horizontally compared between infected and negative control groups at 6, 24, and 48 hpi by DESeq2. A total of 1,897 DEGs were generated between the infected and uninfected samples at 6, 24, and 48 hpi, with 1,260 being upregulated and 637 downregulated ( Table 1).

GO Term Analysis
GO analysis was applied to determine the functions of the DEGs under the thresholds of p < 0.05 and FDR< 0.05.    (Figure 3). Based on the GO terms, genes related to "immune system process, " "response to stimulus" and several other biological processes accounted for comparable percentage in 6 hpi samples (Additional file 7: Table S7). There were 76 DEGs in five GO terms related to "innate immune response" (GO: 0045087), such as immune "system process" (GO: 0002376), "response to stimulus" (GO: 0050896), "immune response" (GO: 0006955), "response to stress" (GO: 0006950), and "defense response" (GO: 0006952) were significantly enriched (p adj < 0.05, Figure 4). A total of 115 DEGs in nine GO terms related to "response to virus" (GO: 0009615) were also significantly enriched in 6 hpi samples infected by SADS-CoV (Figure 4). This indicates that the infection of SADS-CoV could induce the innate immunity and antivirus response in the IPEC-J2 cells. In 24 hpi samples, there were 5,144 genes annotated in 383 GO terms, 4,882 were upregulated, and 262 were downregulated; of the 4,882 upregulated genes, 4,703 (96.33%) were mapped into FIGURE 3 | Functional map of differentially expressed genes enriched for GO terms. All categories were statistically significant (p adj < 0.05). The bars represent the -log 10 p adj -value of the comparison of respective GO terms from infected and uninfected samples. Red, green, and blue bars indicated GO terms clustered in the biological process, cellular component, and molecular function terms, respectively. the BP terms, 111 (2.15%), and 75 (1.52%) were mapped into the CC and MF terms, respectively, (Additional file 8: Table S8). A total of 277 DEGs in seven GO terms were significantly enriched in "defense response" (GO: 0006952), "regulation of response to stimulus" (GO: 0048583), and "regulation of immune system process" (GO: 0002682) in 24 h post SADS-CoV infection (Additional file 9: Figure S1A). Furthermore, 113 DEGs in six GO terms were significantly enriched in protein binding, receptor binding, cytokine activity, cytokine receptor binding, growth factor activity, and growth factor receptor binding, which indicates the cytokine-cytokine receptor interactions were significantly upregulated in virus infected IPEC-J2 cells at 24 hpi (Additional file 9: Figure S1B).
In 48 hpi samples, there were 3,533 genes annotated in 232 GO terms, 3,213 were upregulated and 320 were downregulated; of the 3,213 upregulated genes, 3,081 (95.89%) were mapped into the BP terms, 101 (3.14%) and 31 (0.96%) were mapped into the CC and MF terms, respectively, (Additional file 10: Table S10). The most mapped GO item in all of the three time intervals (6, 24, and 48 hpi) was response to stimulus (GO: 0050896). In biological process category, 224 DEGs were significantly enriched in immune system processes GO terms and inflammatory response GO terms (Additional file 11: Figure S2A). In category of molecular function, 96 DEGs were mapped to protein-and receptor-binding, and cytokine-cytokine receptor binding GO terms (Additional file 11: Figure S2B). This indicates more intensify immune response, inflammation, and cytokines were induced as the infection was more serious in 48 h post SADS-CoV infection.

KEGG Pathway Analysis
The 1,897 DEGs were mapped to KEGG Ontology (KO), and 19, 38, and 52 KO categories were clustered in 6, 24, and 48 hpi samples. Noticeably, NOD-like receptor signaling pathway, Cytokine-cytokine receptor interaction, Tolllike receptor signaling pathway, TNF signaling pathway, and Chemokine signaling pathway were significantly enriched signaling pathways, which indicated the signals were promptly induced and intensified responded during the SADS-CoV infection (Figure 5). In 6 hpi samples, Infectious disease: viral (KEGG ID: ssc05164, ssc05162, ssc05168, ssc05160, and ssc05167), immune system, and signal transduction were the predominant terms (Additional file 12: Table S10). In the 24 hpi samples, 631 DEGs were mapped to categories of Environmental Information Processing, Human Diseases, and Organismal Systems in KO. Same as that of 6 hpi samples, Infectious disease: viral, immune system, and signal transduction were the predominant terms (Additional file 13: Table S11). In 48 hpi samples, 937 DEGs were clustered into 4 KO categories, Cellular Processes, Environmental Information Processing, Human Diseases, and Organismal Systems (Additional file 14: Table S12). Interestingly, the PI3K-Akt and Jak-STAT signaling pathways were significantly upregulated in SADS-CoV infected cells.

Validation of DEGs in RNA-Seq by qRT-PCR
A total of 20 DEGs in the RNA-seq results were further verified by qRT-PCR. Our results shown that the mRNA expression level of the 20 randomly selected genes were in agreement with the expression of RNA-seq (Figure 6, Additional file 15: Table S13), which indicated that the RNA-seq data was reliable. Furthermore, the relative mRNA expression of these genes was represented in the Table S7. Overall, the validation of housekeeping and other selected genes by qRT-PCR demonstrated that the RNA-Seq results were reliable.

The Effects of SAA-3 Protein on SADS-CoV Replication
From the RNA-Seq dataset, SAA-3, an acute phase response protein, was significantly upregulated during the 6, 24, and 48 hpi samples (Additional files 4: Table S4, Additional files 5: Table S5, and Additional files 6: Table S6). To evaluate the effects of SAA3 on SADS-CoV infection, the recombinant IPEC-J2-OE-SAA3 cells with over-expression SAA3 and recombinant IPEC-J2-KD-SAA3 with interference of SAA3 were constructed. Western blot analysis showed that SAA3 was over expressed in IPEC-J2-OE-SAA3 cells, and obviously lower expressed in IPEC-J2-KD-SAA3 cells (Figure 7A). Then, the induction of SAA3 protein and TCID50 of SADS-CoV in IPEC-J2-OE-SAA3, IPEC-J2-KD-SAA3, and normal IPEC-J2 cells were compared. The results showed that SADS-CoV replication was significantly inhibited in the over-expression of SAA3. In contrast, shRNA interference silencing SAA3 gene promoted the expression of SADS-CoV in IPEC-J2 cells (Figure 7B).

DISCUSSION
Viral diarrhea caused by CoVs is still a severe problem in swine. Besides the known etiologies of PEDV and PDCoV, the emergence of SADS-CoV has made the situation more complicated (3,8,9,33). Since SADS-CoV firstly was identified in We systematically analyzed the host defense and biological imperative modulations during infection of SADS-CoV in the annotated genes. We observed that many of the pathways and biological processes in SADS-CoV infected IPEC-J2 cells were significantly upregulated or downregulated. In particular, those genes, which were mainly involved in interactions of signal molecules, signal transduction, cell growth and death system, and immune system-related signaling pathways, were upregulated.
The innate immune response relies on recognition of evolutionarily conserved structures on pathogens, termed pathogen-associated molecular patterns (PAMPs), through a limited number of germ line-encoded pattern recognition receptors (PRRs), including Toll-like receptors (TLR), NODlike receptors (NLR), and RIG-I-like receptors (RLRs) (34). The PAMP recognition of the extracellular RNA by PRRs (such as TLR2, TLR3, and TLR9 receptors) stimulated and signaled to the host the presence of infection and trigger inflammatory cytokines and antivirus responses by activating a multitude of intracellular signaling pathways, including adaptor molecules, kinases, and transcription factors. Then, the activation of gene expression and synthesis of a broad range of molecules, up-regulation of IL-6, CXCL2, IL-1β, IL-8, and other cytokines expression. For example, PEDV-infected porcine intestinal epithelial cells induces the rapid activation of the NF-κB pathway via the TLR family genes, while TLR2-, TLR3-, and TLR9silenced genes can significantly inhibit the expression of NF-κB to promote the replication of PEDV (35). NF-κB is a transcription factor which regulates the expression of many factors that involve in immune system stimulation including a variety of pro-inflammatory cytokines, chemokines, and adhesion molecules (36). In this study, we demonstrated that the SADS-CoV infection significantly upregulated the expression of NFKBIA, NFKBIE, NFKBIB, NFKBID, and NFKBIZ, genes coding for inhibitors of κB (IκB) in IPEC-J2 cells during the three time points (6, 24, and 48 dpi), indicating the SADS-CoV infection could inhibit the expression of NF-κB and, thus, create a favorable milieu for SADS-CoV replication. NF-κB is a transcription factor responsible for modulating the expression of many genes involved in innate immunity, cell proliferation, differentiation, apoptosis, and metastasis. NF-κB interacts with IκB inhibitory proteins to regulate gene expression. In resting cells, NF-κB proteins are predominantly cytoplasmic, associating with members of the inhibitory IκB family. IκB proteins were originally thought to sequester NF-κB in the cytoplasm by masking its nuclear localization sequences (NLSs) (37). So, the upregulated of IκB in SADS-CoV infected cells might lead to the downregulated the activity of NF-κB, that inhibits the IFNs production. The antagonist of innate immune to promote viral survival/replication has been observed in many CoVs infections, including SARS-CoV, IBV, TGEV, PEDV, and PDCoV (38,39). Thus, the results of our study are consistent with the previous studies. A recent study revealed  that SADS-CoV interrupted poly (I:C)-induced phosphorylation and nuclear translocation of IRF3 and NF-κB (26). It might be another signal pathway that SADS-CoV antagonizes in terms of IFN production by upregulated IκB to inhibit the NF-κB pathway; further study is needed to illuminate the molecular mechanism.
SAA is one of the main acute phase proteins that are upregulated rapidly in response to infection, inflammation, and tissue damage in vertebrates (40,41). Moreover, SAA induces a series of inflammatory mediators, including IL-1β, TNF-α, and IL-6 by binding to TLR2 and FPR2 receptors, and resulted in strong biological effects. In many virus infections, SAA is upregulated and plays distinct roles in the virus entry and replication in host cells (42). SAA interacts with hepatitis C virus particles to block virus entry into target cells (43). In present study, SADS-CoV replication was significantly inhibited in the over-expression of SAA3, while enhanced by silencing SAA3 in IPEC-J2 cell. Thus, SAA3 showed inhibition on SADS-CoV replication. The mechanism behind SAA3 suppression of the replication of SADS-CoV is still unknown, however, and further study is needed.
The replication of the virus depends on the host cell internal environments, so many viruses have evolved mechanisms that activate the host cell DNA damage signaling pathway, thereby blocking the cycle to produce the appropriate intracellular environment for proliferation (44)(45)(46). SARS coronavirus (SARS-CoV), murine hepatitis virus (MHV), avian infectious bronchitis virus (IBV), porcine transmissible gastroenteritis virus (TGEV), and PEDV can block the cell cycle (18,21,47,48). According to the results in this study, proteins involved in functional cluster of cell cycle were significantly altered in SADS-CoV-infected cells when compared to the negative control.
In conclusion, this is the first report the interaction between SADS-CoV and IPEC-J2 cells by high-throughput RNA-Seq. Comprehensive functional analysis of host mRNA profiles revealed that SADS-CoV infection induced strong immune responses, including innate immunity, and cytokinecytokine receptor interactions. We also identified the upregulated modulation activating the antiviral defenses of host cells through the elevated expression of immune-related genes and the changing the signaling pathways. The findings of this study could provide important insights into the modulation of host metabolism during SADS-CoV infection and have the potential to improve our understanding the pathogenesis of SADS-CoV.

DATA AVAILABILITY STATEMENT
All sample raw reads deposited at the Short Reads Archive (SRA) database belongs to the National Center for Biotechnology Information (NCBI) and are available under Bioproject ID PRJNA622652.