key: cord-1000610-bb8ph8t8
authors: Peng, Ouyang; Wei, Xiaona; Ashraf, Usama; Hu, Fangyu; Xia, Yongbo; Xu, Qiuping; Hu, Guangli; Xue, Chunyi; Cao, Yongchang; Zhang, Hao
title: Genome-wide transcriptome analysis of porcine epidemic diarrhea virus virulent or avirulent strain-infected porcine small intestinal epithelial cells
date: 2022-01-18
journal: Virol Sin
DOI: 10.1016/j.virs.2022.01.011
sha: 2e42e6ef19515bffb8c16af0cfd8a911b13c0aea
doc_id: 1000610
cord_uid: bb8ph8t8

Porcine epidemic diarrhea virus (PEDV) is the main cause of diarrhea, vomiting, and mortality in pigs, which results in devastating economic loss to the pig industry around the globe. In recent years, the advent of RNA-sequencing technologies has led to delineate host responses at late stages of PEDV infection; however, the comparative analysis of host responses to early-stage infection of virulent and avirulent PEDV strains is currently unknown. Here, using the BGI DNBSEQ RNA-sequencing, we performed global gene expression profiles of pig intestinal epithelial cells infected with virulent (GDS01) or avirulent (HX) PEDV strains for 3, 6, and 12 ​h. It was observed that over half of all significantly dysregulated genes in both infection groups exhibited a down-regulated expression pattern. Functional enrichment analyses indicated that the differentially expressed genes (DEGs) in the GDS01 group were predominantly related to autophagy and apoptosis, whereas the genes showing the differential expression in the HX group were strongly enriched in immune responses/inflammation. Among the DEGs, the functional association of TLR3 and IFIT2 genes with the HX and GDS01 strains replication was experimentally validated by TLR3 inhibition and IFIT2 overexpression systems in cultured cells. TLR3 expression was found to inhibit HX strain, but not GDS01 strain, replication by enhancing the IFIT2 expression in infected cells. In conclusion, our study highlights similarities and differences in gene expression patterns and cellular processes/pathways altered at the early-stage infection of PEDV virulent and avirulent strains. These findings may provide a foundation for establishing novel therapies to control PEDV infection.

Coronaviruses (CoVs) are enveloped, positive-sense, single-stranded, RNA viruses that are known to infect a variety of mammalian species, such as humans, pigs, bats, and camels (Haagmans et al., 2014; Niederwerder et al., 2018; Wong et al., 2007; Zhou et al., 2020) . CoVs, belonging to the family Coronaviridae, are classified into four genera: Alphacoronavirus, Betacoronavirus, Gammacoronavirus, and Deltacoronavirus Woo et al., 2012) . Of these, the viruses belonging to the genus Betacoronavirus can jump from animals to humans, and thus, can be responsible for lethal acute respiratory illnesses in humans such as the Middle East respiratory syndrome and the severe acute respiratory syndrome (Drosten et al., 2003; Zaki et al., 2012) .

Notably, the recent pandemic of Coronavirus Disease 2019, caused by the severe acute respiratory syndrome-CoV2, infected over 80 million people with exceeding 1.8 million deaths globally (Sun et al., 2020; WHO, 2021; Wu et al., 2020) . Moreover, the bat-HKU2-like novel porcine CoV has been identified to cause acute diarrhea syndrome in pigs (Gong et al., 2017; Zhou et al., 2018) . The replication of bat-HKU2-like CoV in multiple mammalian cells highlighted its ability to cross the species barrier (Edwards et al., 2020) , which requires great attention.

Porcine epidemic diarrhea virus (PEDV), belonging to the genus Alphacoronavirus, is the causative agent of porcine epidemic diarrhea (PED) disease, and it has been a considerable threat to pigs of all ages with a mortality rate of~100% in suckling piglets caused by virulent strains (Li et al., 2012) . PED is a highly devastating enteric disease which is characterized by vomiting, diarrhea, dehydration, weight loss, and eventually death . Since its first isolation in the United Kingdom and Belgium in the 1970s, PEDV has been reported in several countries, including the United States, Canada, Mexico, China, Korea, and Japan, resulting in serious damage to the pig industry and huge economic losses (Chen et al., 2012; Chung et al., 2015; Mole, 2013; Niederwerder et al., 2018; Takahashi et al., 1983) . PEDV primarily infects the cells of the small intestine, and pig-to-pig transmission mainly occurs by fecal-oral route . PED outbreaks were found to be associated with the novel virulent recombinant PEDV strains (G2 subgroup), which were genetically distinct from the classic avirulent CV777 strain (G1 subgroup) Huang et al., 2013; Wang et al., 2020) .

High-throughput transcriptomic approaches are extensively used for studying the molecular forces deriving the virus-host systems. The RNAsequencing analysis of the PEDV-infected cultured pig intestinal cells revealed the perturbation of signal transducer and activator of transcription, cell cycle, and apoptosis pathways at the later stages of infection . In a similar context, a few studies have identified several key host genes involved in regulating the antiviral immune response against PEDV such as FOSL1, IL1A, ISG15, and some other ISGs Sun et al., 2019; Wang et al., 2019) . However, the host transcriptomic landscape at an early phase of PEDV infection and differences in the host response upon infection with an avirulent or a virulent PEDV strain are currently unknown.

In this study, using the DNBSEQ RNA-sequencing, we comprehensively analyzed the genome-wide transcriptomic profiles of cellular mRNAs isolated from cultured porcine small intestinal epithelial cells (IPEC-J2 cells) infected with virulent (GDS01) or avirulent (HX) PEDV strains at three points of early stages of virus infection. Our data identified PEDV-induced widespread alterations in cellular mRNA expression and highlighted key differences in the host response to a virulent and an avirulent PEDV strain. This study may shed light on the biology of PEDV infection that can be considered for developing antiviral therapies against PEDV.

Porcine intestinal columnar epithelial cells (IPEC-J2 cell line) were grown in Dulbecco's modified Eagle's medium (#11995, Gibco, New York NY, USA) containing 5% fetal bovine serum (#10100, Gibco), 100 U/mL penicillin, and 100 mg/mL streptomycin. PEDV strains HX (avirulent, accession number: MH726408) and GDS01 (virulent, accession number: KM089829), belonging to the subgroup G1 and G2, respectively, were isolated and propagated in our laboratory, which were described by our previous study . The median tissue culture infectious dose (TCID 50 ) of both viral strains were as follows: HX strain, 1.47 Â 10 À5 /0.1 mL and GDS01 strain, 2.15 Â 10 À6 /0.1 mL.

IPEC-J2 cells plated in six-well plates (1 Â 10 6 cells/well) were infected with PEDV strain GDS01 or HX at an MOI of 1 for a period of 0, 3, 6, and 12 h. After infection, cells were washed with 1 Â phosphate- A Workflow. IPEC-J2 cells were mockinfected or infected with PEDV strain GDS01 or HX (MOI ¼ 1), followed by sample collection at 3, 6, 12 hpi. Samples from each group were prepared in triplicate. Total RNA was extracted and BGI DNBSEQ RNAsequencing was performed. B Immunofluorescence staining of the viral S-protein in PEDV-infected IPEC-J2 cells. Cells were fixed at 3, 6, or 12 hpi, and S protein (green) was detected by indirect immunofluorescence assay. Nuclei (blue) were shown by 4 0 ,6-diamidino-2-phenylindole (DAPI) staining. The images of cells were acquired by a fluorescence microscope (Nikon Eclipse 80i) at a 20Â magnification. C Viral titers in the culture supernatants were determined by TCID 50 assay. Data are shown as mean AE SEM and representative of three independent experiments. Data were analyzed by the Mann-Whitney test. *P < 0.05. PEDV, porcine epidemic diarrhea virus; MOI, multiplicity of infection; SEM, standard error of mean. buffered saline (PBS) followed by total RNA extraction using the TRIzol reagent (#15596, Thermo Fisher Scientific, Waltham MA, USA). Samples harvested at 0 h post-infection (hpi) were considered as mock-infected. Three biological replicates were set per time point (Fig. 1A ).

The mRNA of samples was purified by magnetic beads connected with oligo(dT) and were then fragmented into small pieces with the fragment buffer. Random hexamer-induced reverse transcription was used to generate first-strand cDNA, followed by second-strand cDNA synthesis. Afterward, the A-tailing mixture and RNA Index Adapters were added by incubation to end repair. The cDNA fragments obtained from the previous step were amplified by PCR, and products were purified by Ampure XP beads, and then dissolved into the EB solution. The products were evaluated using the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara CA, USA). The double-stranded PCR products from the previous step were heated, denatured, and circularized by the splint oligo sequence to obtain the final library. The final library was then amplified with phi29 to prepare DNA nanoball (DNB) with more than 300 copies of one molecule. DNBs were loaded into the patterned nanoarray and single-ended 150 base reads were generated on the BGIseq500 platform (BGI, Wuhan, China).

The quality control of raw reads was performed using the SOAPnuke software (Cock et al., 2010) . The following reads were filtered out: 1) the reads containing adaptors, 2) the reads with N content more than 5%, and 3) low-quality reads (reads with a base quality score of less than 10 accounted for more than 20% of the total bases as low-quality reads). The filtered clean reads obtained from each sample were~45.83 million, packaged in the form of FASTQ, with an average size of 6.87 Gb. The Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT) was used to map RNA-sequencing reads to the porcine reference genome (RefSeq assembly accession: GCF_000003025.6). Firstly, we use HISAT global FM index to anchor the position of partial sequences in each read on the genome, and then, use the partial genome indexes of these alignment positions to align the remaining sequences of each read to extend the alignment area .

In order to quantitatively analyze the gene expression level, Bowtie2 (V2.2.5) was employed to align the clean reads to the porcine reference gene sequences (Langmead et al., 2012) . The soft parameters were set as: -q -phred64 -sensitive -dpad 0 -gbar 99999999 -mp 1,1 -np 1 -score-min L, 0, À0.1 -p 16 -k 200. And then, the RSEM (V1.2.8) was used to calculate and normalize the gene expression level in each sample with the default parameters . Differentially expressed genes (DEGs) were selected by applying the following filtering parameters: | log2 (fold change)| ! 1, and Q-value 0.01.

GO term and pathways enrichment analysis were performed respectively using GO::TermFinder (Boyle et al., 2004) and Kyoto encyclopedia of genes and genomes (KEGG) pathways database (Kanehisa et al., 2008) . For GO term analysis, we mapped all DEGs to each term in the GO database (Harris et al., 2004) , calculated the number of genes in each term, and then applied the hypergeometric test to find GO terms that were significantly enriched in candidate genes compared with the background of all genes of the pig. The calculated P-value was corrected by Bonferroni. The calculation method for predicting enriched pathways was the same as that described for the GO terms. The GO term or pathway with Q-value (corrected P-value) 0.05 was defined as a significantly enriched GO term or pathway corresponding to the candidate gene.

To analyze mRNA expression, 0.1 μg of total RNA was used to reverse transcribe in order to construct cDNA using the First Strand cDNA Synthesis Kit (#FSK-101, TOYOBO, Tokyo, Japan) containing oligo(dT) primer (5 0 -TTTTTTTTTTTTTTTTTTTT-3 0 ) according the manufacturer's protocol. Quantitative real-time PCR was performed using the SYBR Green Real-time PCR Master Mix (#11201ES03, Yeasen, Shanghai, China) and a LightCycler480 II system (Roche, Basel, Switzerland). Amplification was performed at 50 C for 2 min, 95 C for 10 min, and then followed by 40 cycles of 95 C for 15 s, 60 C for 15 s, and 72 C for 30 s. The relative expression values of candidate mRNAs were normalized to that of GAPDH in each sample using the 2 ÀΔΔCt method.

The number of PEDV genomic copies was assessed with a SYBR Green Real-time PCR Master Mix (Yeasen) and a LightCycler480 II system (Roche). The qRT-PCR conditions used were a holding stage of two steps, which included a first step at 48 C for 30 min, and a second at 95 C for 10 min. The cycling stage included 40 repeated cycles of two steps, a first step of 15 s at 95 C and a second step of 1 min at 60 C. The melt curve stage was composed of four steps, a first step of 1 min at 95 C, a second step of 1 min at 60 C, a third step with a gradual increase in temperature with 0.35 C for 0.3 s to obtain a temperature of 95 C and a fourth step of 15 s at 60 C. The primers used for the PCR were designed from the conserved regions of the PEDV nucleocapsid gene for universal detection of the strain used for inoculation (forward, 5 0 -CGCAAA-GACTGAACCCACTAA-3 0 ; reverse, 5 0 -TTGCCTCTGTTGTTACTTGGAGAT-3 0 ).

IPEC cells were seeded in 6-well plates at a density of 1 Â 10 6 cells/ well and grew for 24 h. Non-adherent cells and media were removed and replaced with fresh DMEM, and adherent IPEC cells were then treated with commercially available selective TLR3 inhibitor CU CPT 4a (#4883, R&D, Minneapolis MN, USA) at IC 50 ¼ 3.44 μM for a period of 1 h. Later, treated cells were mock-infected or infected with GDS01 or HX strain (MOI ¼ 1) for a period of 3, 6, and 12 h. Subsequently, samples were collected and processed for further experiments.

For the construction of the pCMV-Flag-IFIT2 plasmid, the cDNA derived from IPEC cells was utilized as a template to amplify the porcine IFIT2 gene by PCR, followed by cloning into the pCMV-Flag vector ( Supplementary Fig. S3 ). IPEC cells were seeded in 6-well plates at a density of 1 Â 10 6 cells/well and grew for 24 h. Cultured cells were either transfected with empty vector or pCMV-Flag-IFIT2 (2 μg each) using the Lipofectamine 3000 (#L3000008, Thermo Fisher Scientific). At 36 h post-transfection, the cells were infected with GDS01 or HX strain (MOI ¼ 1) for a period of 3, 6, and 12 h. Subsequently, samples were collected and processed for further experiments.

IPEC cells were seeded in 6-well plates as described above and incubated for 24 h. Then the culture medium was removed and polyinosinicpolycytidylic acid [Poly(I:C)] was transfected at 0.5 mg/mL using Lipofectamine 3000. At 12 h post-transfection, cells were infected with GDS01 or HX strain at an MOI of 1 for 3, 6, and 12 h. Subsequently, the cell supernatants were collected and processed for the quantification of viral copies and titers by quantitative real-time PCR and TCID 50 assay, respectively.

To investigate the dynamics of host genes expression during the early stages of the PEDV infection, cultured IPEC-J2 cells were mock-infected or infected with either the PEDV GDS01 strain (virulent) or PEDV HX strain (avirulent) at MOI of 1 for a period of 3, 6, and 12 h (Fig. 1A) . To validate whether IPEC-J2 cells have acquired successful infection, the replication of GDS01 and HX viruses was examined by an immunofluorescence assay (IFA) and the TCID 50 assay. It was observed that both viral strains replicated successfully in a time-dependent manner as assessed by the increasing number of PEDV S protein-positive cells (Fig. 1B) and the increasing production of infectious viral particles (Fig. 1C) . In both assays, the GDS01 strain was observed to reach at significantly higher virus titers compared to the HX strain at 6 h and 12 h infection time-points ( Fig. 1B and C) . Overall, these results indicate the establishment of PEDV GDS01 and HX strains infection in cultured IPEC-J2 cells.

To examine the reproducibility and the specificity of each group, the gene expression levels were used to conduct a principal component analysis (PCA) for each biological replicate. Every sample from the same group was clustered closely together, which suggested that the reproducibility of each treatment was satisfactory, and the specificity between groups was apparent ( Fig. 2A) .

In order to visualize the transcriptomic profile of the GDS01-or HX -infected cells, a total of 21 sequencing libraries prepared at 0, 3, 6, and 12 hpi were sequenced, and~976 million raw reads were obtained. The number of clean reads ranged from~44 to 48 million after filtering out the adapter and low-quality reads. To this end, a total of 1965 significant DEGs were identified in infected samples when compared to mock infected samples under the parameter of Q-value 0.01 and |log2 (fold change)| ! 1. Overall, approximately half of the genes in HX-infected groups at 3 hpi are down regulated, whereas genes in other groups exhibiting the down-regulated expression pattern constitute two third of total DEGs identified, and the increase in the number of down-regulated genes is dependent on the time of infection (Fig. 2B) . Moreover, the total number of DEGs detected in all three infection time-points is more in the GDS01-infected group compared to the HX-infected group (Fig. 2B) , which may associate with the differences in the replication of both viral strains as observed in Fig. 1B and C. Taken together, these data demonstrate that the infection of cultured IPEC-J2 cells by PEDV GDS01 or HX strain induces widespread alterations in the expression pattern of host genes.

Venn diagrams were generated to examine the unique and overlapping DEGs among subgroups with the same virus strain infection at three infection time-points or among subgroups with different strains infection but the same infection time-point. To this end, 55 and 38 DEGs were noticed common among all three subgroups in the GDS01 and HX groups, respectively. However, no DEG was found common among all three infection time-points belong to both infection groups, suggesting a marked differential host response at each stage of early infection (Supplementary Fig. S1B ). In the HX subgroups, immune-associated DEGs, such as GATA3, FOS, IFIR2, IL1R2, CXCL10, were detected at each timepoint, whereas in the GDS01 subgroups, only IL1R2 and CXCL10 genes were identified.

The expression values of 12,604 genes detected in all 21 samples were used to construct the co-expression module using a gene hierarchical cluster network tool, R package, WGCNA (V1.69) (Peter et al., 2008) . The correlation matrix and the adjacency matrix of the gene expression profile of the PEDV (GDS01 or HX)-infected samples were analyzed through the standard parameters of WGCNA. The co-expression modules were constructed and the hierarchical average linkage clustering method was utilized to identify the gene modules of each gene network. All samples were clustered, and a total of 16 different module eigengenes were observed based on their gene expression level (Fig. 3A) . The total number of genes in all 16 modules is shown in Table 1 . Interactions of the 16 co-expression modules were also analyzed and were depicted as a heatmap (Fig. 3B) .

To determine how IPEC-J2 cells respond to the PEDV infection, 16,942 and 16,000 genes identified in the HX and the GDS01 infection group, respectively, were classified into a total 12 of clusters (six clusters to each infection group) using the Mfuzz package (V2.48.0) (Lokesh et al., 2007) by considering their expression tendencies along with the Fig. 4A and B) . In the GDS01 group, the genes classified in clusters 1, 2, and 5 were down-regulated, whereas the genes in clusters 3, 4, and 6 demonstrated up-regulation. Similarly, in the HX groups, clusters 1, 4, and 6 presented up-regulated genes, and clusters 2, 3, and 5 showed down-regulated genes. Those genes that clustered together were more likely considered to be classified into the same functional gene set.

To know the biological importance of the DEGs, GO term enrichment analysis was performed by considering three functional categories that included molecular function, cellular components, and biological processes. For the stringency purpose, the top 30 or less highly enriched GO terms (Q-value < 0.05) corresponding to all DEGs were considered (Fig. 5) . At each infection time-point, the GDS01 and HX groups shared a similar GO enrichment pattern, but the number of DEGs enriched in GDS01 group was more than that in HX group. The genes enriched in the molecular function category are mainly associated with DNA-binding and DNA-binding transcription factor activity. The genes enriched in cellular component were markedly less than the other two GO categories, and were predicted to localize primarily in the nucleus, cell periphery, and intermediate filament. The genes relevant to nuclear distribution were enriched at 3 and 12 hpi in the GDS01 group but only at 12 hpi in the HX group. The biological process classification predicted the enrichment of DEGs mainly in the regulation of cellular transcription, RNA polymerase II activity, cellular proliferation, and apoptosis.

To identify the cellular pathways potentially involved during the PEDV infection, the KEGG database was employed. The top 10 highly enriched biological pathways corresponding to each condition in both Fig. 3 . Construction of a weighted gene co-expression network. A A cluster dendrogram was built based on the dis-similarity of the topological overlap, which presented 16 gene co-expression modules. The grey module indicates no co-expression between the genes. B Correlated heatmap plot of the adjacency modules in the WGCNA. The rectangle of each row and column represents a module eigengene. In the correlated heatmap plot, light blue represents low adjacency, whereas red represents high adjacency. WGCNA, weighted correlation network analysis. infection groups are shown together in a bubble diagram (Fig. 6 ). Upon infection with the HX strain, the KEGG pathways were observed to be enriched in the immune-associated pathways such as TNF and IL-17 signaling. In contrast, the GDS01 infection resulted in marked induction of apoptosis-and autophagy-related pathways, including PI3K-Akt, p53, and FoxO signaling. The most commonly enriched pathways among all the infection conditions are Hippo signaling and osteoclast differentiation. In addition, the highest number of DEGs were found to associate with the PI3K-Akt signaling in the GDS01 group and with the MAPK signaling in the HX group. In order to understand the host immune defense response upon PEDV infection, the four most enriched immune signaling pathways are (Fig. 7B) . The upregulation of TLR3, PI3K and the down-regulation of TRAIL, and RIP1 were noticed in a time-dependent manner in both infection groups. Two DEGs exhibited a time-specific expression pattern in both infection groups: reduction of FADD at 12 hpi and elevation of JunB at 6 hpi. Moreover, LCN2/PEPCK/TNFR2 and LIF/FOS showed almost constant up-and down-regulated expression, respectively, across all three infection time-points studied in both groups.

Overall, these findings predict that both PEDV strains GDS01 and HX may induce the alteration of the host cellular pathways to facilitate their own replication.

To confirm the differential expression of the PEDV-induced genes detected in our sequencing data, we chose four DEGs (TLR3, TRAIL, PI3K and JunB) as predicted above to be engaged in immune response, apoptosis, and autophagy, and then validated their expression in PEDV (GDS01 or HX)-infected or mock-infected IPEC-J2 cells by quantitative real-time PCR. Primers used in our study are listed in Table 2 . As shown in Fig. 8 , these DEGs exhibited similar expression signatures as were detected in the sequencing data, indicating the robustness of our experimental settings and bioinformatic analyses. When compared to the gene expression pattern upon HX infection, TLR3 was observed to show an upregulation pattern at three infection time-point, TRAIL and TLR6 genes exhibited a downregulated pattern at 12 h infection time-point, whereas PI3K and JunB showed an upregulated pattern at 6 h infection timepoint, upon GDS01 infection. Overall, these results validate our sequencing data and highlight some viral strain-dependent gene expression patterns.

To discern a possible role of detected genes in regulating the PEDV replication, we selected TLR3 as one of the DEGs in the sequencing data. In the HX-infected IPEC-J2 cells, the TLR3 expression was observed to be significantly increased than that in the GDS01-infected cells, which might have contributed to the lower propagation of HX strain when compared to GDS01 strain (Fig. 7B) . To verify this speculation, we employed a selective TLR3 inhibitor in a series of experiments and examined its effect on the downstream antiviral IFIT2 gene and PEDV (HX and GDS01) replication in IPEC-J2 cultured cells. It was observed that TLR3 inhibitor caused a reduction of IFIT2 mRNA expression level (Fig. 9A) . Furthermore, the reduced expression of IFIT2 mRNA upon TLR3 inhibition was associated with significantly enhanced viral copies and titers of HX strain, as measured by the quantitative real-time PCR and TCID 50 assay, respectively ( Fig. 9B and C) . No significant effect of TLR3 inhibition was noticed on the viral copies and titers of GDS01 strain ( Fig. 9B and C) . We then examined the effect of poly(I:C) treatment on the replication of both viral strains. It was observed that the replication of HX strain in IPEC cells was significantly suppressed upon poly(I:C) stimulation, whereas the replication of GDS01 strain was impervious to poly(I:C) transfection, as determined by the quantitative real-time PCR and TCID 50 assay ( Fig. 9D and E) . Hence, these data suggest a GDS01specific immune escape mechanism that facilitates the GDS01 strain replication.

We next confirmed the role of the IFIT2 genes in the TLR3-mediated inhibition of HX strain. To this end, cultured IPEC-J2 cells were transfected with empty vector or Flag-IFIT2 plasmids, followed by infection with HX or GDS01 strain for the indicated time points. The efficient transcription of IFIT2 mRNAs in transfected cells was validated by quantitative real-time PCR (Fig. 10A) . Overexpression of IFIT2 caused a marked suppression of viral copies and titers of both HX and GDS01 strains, determined by quantitative real-time PCR and TCID 50 assay, respectively ( Fig. 10A and B) .

Taken together, these data indicate that TLR3, via augmenting the antiviral IFIT2 gene expression, restricts the HX strain replication.

Advancements in RNA-sequencing technologies have provided us an opportunity to understand the complex host-virus interaction in detail (Jones et al., 2014; Liu et al., 2019) . Previous studies have performed the global analysis of host transcriptomic signatures to gain insights into the host response to PEDV infection at the late stages of virus infection Shen et al., 2020) . Herein, we provide, for the first time, a view of genome-wide alterations in the host transcriptome that occur in response to PEDV infection, with a focus on differences in the host response to virulent (GDS01) and avirulent (HX) PEDV strains at the early stages of infection. Upon RNA-sequencing analysis of cultured IPEC-J2 cells infected with either GDS01 or HX viruses, we detected a total of 1968 significant DEGs in both infection groups (GDS01 group: 1065 genes and HX group: 903 genes), which were predicted to regulate mainly cellular transcription, antiviral immune response, autophagy, and apoptosis. The findings obtained from this study highlighted the intricate regulation of cellular genes during the initial phases of PEDV infection, and suggested their potential role in regulating the PEDV pathogenesis. Further studies are required to understand the biological impact of identified DEGs during PEDV infection. Whether or not, targeting these genes confers a therapeutic effect against PEDV infection, remains to be evaluated in future studies.

The differences in the gene expression pattern observed upon infection with virulent and avirulent PEDV strains could pave a way to better understand the strain-specific host responses, which could be considered for developing the novel antiviral strategies. Previous studies mainly examined the transcriptomic profile of pig enterocytes infected with a classical avirulent PEDV CV777 strain or Vero E6 cells infected with mutated virulent JS201603 strain (Zhang et al., 2018a) . In contrast, we examined and compared the transcriptomic landscape of pig The comparative analysis of host response to both viral strains revealed that the overall modulation of genes enriched in immunity/inflammation was stronger upon infection with avirulent HX strain, whereas the virulent GDS01 strain intensified the response of autophagy and apoptosis-associated genes, thus uncovering the potential mechanisms of virulence associated with the GDS01 strain. In a recent study, Hu et al. found that the DEGs detected in CV777-infected IPEC-J2 cells were enriched to the inflammatory response (GO:0006954), immune response (GO:0006955), regulation of inflammatory response (GO:0050727), innate immune response (GO:0045087) categories , which was in agreement with our GO term enrichment analysis of DEGs in the avirulent HX group. It would be intriguing to perform the extensive functional analysis of these detected genes in orchestrating the immunity to PEDV infection. Autophagy has been shown to facilitate the replication of a virulent PEDV CH/YNKM-8/2013 stain at the late stages of infection but not at the early stages, via engaging the autophagy regulators and RNA interference in Vero cells (Guo et al., 2017) . The quantitative proteomic analysis of PEDV CV777 vaccine strain-infected Vero cells at 48 hpi indicated that the differentially expressed proteins were highly enriched in the autophagy pathway (Sun et al., 2015) . These results are inconsistent with our finding of strong enrichment of genes in the autophagy pathway at the early stages of virulent GDS01 strain infection in IPEC-J2 cells. Such observations might be associated with the infection stage, strain, and cell-type specificity, which required additional investigations. Furthermore, the proteomic analysis is required to confirm, whether or not, the dysregulated proteins during the early stage of GDS01 strain infection are enriched in the autophagy pathway.

PEDV has been shown to induce apoptosis in vitro and in vivo, which plays an important role in cytotoxicity and pathogenesis (Kim et al., 2014; Xu et al., 2019; Zeng et al., 2015) . For instance, the PEDV virulent SM98-1 strain infection promoted the nuclear translocation of mitochondrial apoptosis-inducing factor, which subsequently enhanced the virus replication in cultured Vero cells and enterocytes of experimentally-infected piglets by inducing apoptosis in a caspase-independent manner (Oh et al., 2020) . In a similar context, the S1 spike protein of virulent SM98-1 strain caused apoptosis of Vero cells, which was found to be associated with the cleavage of apoptosis-inducing factor mitochondria associated-1 and poly (ADP ribose) polymerase-1 and with the activation of caspase-3 and caspase-8 . These findings support our observation of strong enrichment of DEGs in apoptotic pathway during infection of IPEC-J2 cells with the virulent GDS01 strain. Since little is known about the mechanisms by which PEDV induces apoptosis, the genes and pathways identified to be involved in the apoptosis in our study, could be utilized to uncover the precise mechanisms of PEDV-induced apoptosis.

PEDV evades the host immunity by two distinct mechanisms: 1) protecting the viral genome from the cellular nucleic acid sensors and 2) producing the virus-encoding interferon antagonists . PEDV nonstructural protein 1, 3, 5, 7, 8, 14, and 15 blocks the type I and/or type III interferon production by inhibiting the nuclear translocation of NF-κB and IRF1, cleaving NEMO, deubiquitinating RIG-I and STING, or through some other unknown mechanisms Zhang et al., 2017 Zhang et al., , 2018b . PEDV AJ1102 strain nucleocapsid protein sequesters the IRF3-TBK1 interaction, which leads to the suppressed activity of NF-κB and the reduced synthesis of type I interferon (Ding et al., 2014) . The difference in the activation or suppression of immune signaling pathways has also been associated with the strain-type, time of infection, or cell-specificity (Du et al., 2019) . For example, PEDV virulent non-S-INDEL strain impedes the NF-κB signaling pathway and reduces subsequent cytokine/interferon production by negatively regulating the TLR4, TLR7, TLR8, and TLR9 expression, whereas the less-virulent PEDV S-INDEL strain enhances the cytokines production through the non-canonical NF-κB pathway by stimulating RIG-I (Temeeyasen et al., 2018) . In our study, we observed that the HX strain stimulated the TLR3 and the downstream antiviral IFIT2 gene expression to a relatively higher extent than GDS01 strain, which led to reduced HX strain, but not GDS01 strain, replication, highlighting an immune evasion mechanism that was specific to GDS01 strain. Additionally, we observed that the ectopic expression of IFIT2 in cells suppressed both HX and GDS01 replication. The reduced replication of GDS01 strain upon expression of IFIT2, but not TLR3, suggested the engagement of TLR3-independent antiviral mechanisms, which could be intriguing to examine in future studies. Moreover, the obvious differences in the secondary and tertiary RNA structures of both viral strains may contribute to a differential TLR3 recognition, thus causing a strain-specific response.

In summary, we performed the global gene expression profiles of pig enterocytes infected by virulent or avirulent PEDV strains. Upon integration of RNA-sequencing, functional enrichment analyses, gene coexpression network analysis, experimental validations, and functional studies, we highlighted similarities and differences in gene expression pattern and cellular processes/pathways perturbed at the early-stage Fig. 10 . TLR3 induced antiviral IFIT2 gene expression restricts the HX strain replication. A IFIT2 mRNA expression levels measured by quantitative real-time PCR. B Viral copies in the culture supernatants detected through the quantitative real-time PCR. C Viral titers in the culture supernatants determined by TCID 50 assay. Data are expressed as mean AE SEM from three independent experiments. Data were analyzed using the Mann-Whitney test. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. PCR, polymerase chain reaction; TCID 50 , the 50% tissue culture infectious dose; SEM, standard error of mean.

infection of PEDV virulent and avirulent strains. This study may shed light on molecular mechanisms involved during the PEDV replication. As the transcriptomic signatures at the early stages of PEDV infection have not yet been explored, our study may serve as a valuable resource for PEDV-related research studies in the future.

The raw data of RNA-sequencing have been uploaded to National Center for Biotechnology Information (NCBI) (https://www.ncbi.nl m.nih.gov) and the accession numbers are from SRR14552484 to SRR14552504.

This article does not contain any studies with human or animal subjects performed by any of the authors.

Ouyang Peng: conceptualization, software, formal analysis, resources, data curation, writing-original draft preparation, visualization. Xiaona Wei: conceptualization, methodology, investigation, resources. Usama Ashraf: methodology, writing-review and editing. Fangyu Hu: formal analysis, resources, data curation. Yongbo Xia: investigation. Qiuping Xu: writing-review and editing. Guangli Hu: investigation. Chunyi Xue: validation, writing-review and editing, supervision. Yongchang Cao: validation, supervision, funding acquisition. Hao Zhang: conceptualization, methodology, formal analysis, resources, project administration, funding acquisition.

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Complete genome sequence of a porcine epidemic diarrhea virus variant

GO::TermFinder-open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes

Isolation and characterization of porcine epidemic diarrhea viruses associated with the 2013 disease outbreak among swine in the United States

Porcine epidemic diarrhea virus s1 protein is the critical inducer of apoptosis

Isolation of porcine epidemic diarrhea virus during outbreaks in South Korea

The sanger fastq file format for sequences with quality scores, and the solexa/illumina fastq variants

Origin and evolution of pathogenic coronaviruses

Porcine epidemic diarrhea virus nucleocapsid protein antagonizes beta interferon production by sequestering the interaction between IRF3 and TBK1

Identification of a novel coronavirus in patients with severe acute respiratory syndrome

Manipulation of intestinal antiviral innate immunity and immune evasion strategies of porcine epidemic diarrhea virus

Swine acute diarrhea syndrome coronavirus replication in primary human cells reveals potential susceptibility to infection

A new bat-hku2-like coronavirus in swine

Porcine epidemic diarrhea virus induces autophagy to benefit its replication. Viruses-Basel 9

Middle east respiratory syndrome coronavirus in dromedary camels: an outbreak investigation

The Gene Ontology (GO) database and informatics resource

Transcriptome analysis reveals modulation of the stat family in pedv-infected ipec-j2 cells

Origin, evolution, and genotyping of emergent porcine epidemic diarrhea virus strains in the United States

RNA-seq analysis of host and viral gene expression highlights interaction between varicella zoster virus and keratinocyte differentiation

Kegg for linking genomes to life and the environment

Porcine epidemic diarrhea virus induces caspase-independent apoptosis through activation of mitochondrial apoptosis-inducing factor

Hisat: a fast spliced aligner with low memory requirements

Fast gapped-read alignment with bowtie 2

RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome

New variants of porcine epidemic diarrhea virus

An alternative pathway of enteric pedv dissemination from nasal cavity to intestinal mucosa in swine

Porcine epidemic diarrhea virus and the host innate immune response

Rna-seq based transcriptome analysis during bovine viral diarrhoea virus (BVDV) infection

Mfuzz: a software package for soft clustering of microarray data

Pathogenesis of porcine epidemic diarrhea virus isolate (us/Iowa/18984/ 2013) in 3-week-old weaned pigs

Deadly pig virus slips through us borders

Swine enteric coronavirus disease: a review of 4years with porcine epidemic diarrhoea virus and porcine deltacoronavirus in the United States and Canada

Caspase-mediated cleavage of nucleocapsid protein of a protease-independent porcine epidemic diarrhea virus strain

WGCNA: an R package for weighted correlation network analysis

Porcine epidemic diarrhea virus infection blocks cell cycle and induces apoptosis in pig intestinal epithelial cells

Analysis of protein expression changes of the vero e6 cells infected with classic pedv strain cv777 by using quantitative proteomic technique

Transcriptomic analysis of small intestinal mucosa from porcine epidemic diarrhea virus infected piglets

Covid-19: epidemiology, evolution, and crossdisciplinary perspectives

An outbreak of swine diarrhea of a new-type associated with coronavirus-like particles in Japan

Differential gene modulation of pattern-recognition receptor tlr and rig-i-like and downstream mediators on intestinal mucosa of pigs infected with pedv non s-indel and pedv s-indel strains

Global mapping of H3K4 trimethylation (H3K4me3) and transcriptome analysis reveal genes involved in the response to epidemic diarrhea virus infections in pigs

Pathogenicity and immunogenicity of a new strain of porcine epidemic diarrhea virus containing a novel deletion in the n gene

PEDV enters cells through clathrin-, caveolae-, and lipid raft-mediated endocytosis and traffics via the endo-/ lysosome pathway

Bats as a continuing source of emerging infections in humans

Discovery of seven novel mammalian and avian coronaviruses in the genus deltacoronavirus supports bat coronaviruses as the gene source of alphacoronavirus and betacoronavirus and avian coronaviruses as the gene source of gammacoronavirus and deltacoronavirus

A new coronavirus associated with human respiratory disease in China

Porcine epidemic diarrhea virus infections induce apoptosis in vero cells via a reactive oxygen species (ros)/p53, but not p38 mapk and sapk/jnk signalling pathways

Isolation of a novel coronavirus from a man with pneumonia in Saudi Arabia

Proteome analysis of porcine epidemic diarrhea virus (PEDV)-infected vero cells

Inhibition of nf-kappa b activity by the porcine epidemic diarrhea virus nonstructural protein 1 for innate immune evasion

Genome-wide analysis of differentially expressed genes and the modulation of pedv infection in vero e6 cells

Type III interferon restriction by porcine epidemic diarrhea virus and the role of viral protein nsp1 in IRF1 signaling