Differing coronavirus genres alter shared host signaling pathways upon viral infection

Coronaviruses are important viral pathogens across a range of animal species including humans. They have a high potential for cross-species transmission as evidenced by the emergence of COVID-19 and may be the origin of future pandemics. There is therefore an urgent need to study coronaviruses in depth and to identify new therapeutic targets. This study shows that distant coronaviruses such as Alpha-, Beta-, and Deltacoronaviruses can share common host immune associated pathways and genes. Differentially expressed genes (DEGs) in the transcription profile of epithelial cell lines infected with swine acute diarrhea syndrome, severe acute respiratory syndrome coronavirus 2, or porcine deltacoronavirus, showed that DEGs within 10 common immune associated pathways were upregulated upon infection. Twenty Three pathways and 21 DEGs across 10 immune response associated pathways were shared by these viruses. These 21 DEGs can serve as focused targets for therapeutics against newly emerging coronaviruses. We were able to show that even though there is a positive correlation between PDCoV and SARS-CoV-2 infections, these viruses could be using different strategies for efficient replication in their cells from their natural hosts. To the best of our knowledge, this is the first report of comparative host transcriptome analysis across distant coronavirus genres.

Presently, the origins of SARS-CoV-2 are unclear. However, the 2002 SARS-CoV-1 outbreak-and the entry of SARS-CoV-2 into the human population-has been implicated on the wild animal handling practices common in Southern China 3 .
Gamma-and delta-coronaviruses primarily infect birds, but cases of mammalian infections have been documented 1,[12][13][14][15] . CoVs can also have devastating effects in livestock populations, particularly pigs. CoVs with swine health implications include transmissible gastroenteritis virus (TGEV), porcine epidemic diarrhea virus (PEDV), porcine deltacoronavirus (PDCoV), and swine acute diarrhea syndrome virus (SADS-CoV), among others 1,16,17 . The Deltacoronavirus genus comprises mostly avian CoV pathogens of songbirds including HKU11 (bulbul coronavirus), HKU12 (thrush coronavirus), and HKU13 (munia coronavirus) 18 . PDCoV is an emerging viral disease of swine circulating globally with mortality in up to 40% of infected neonatal pigs 19 . PDCoV's usage of an interspecies conserved amino acid domain within aminopeptidase N (APN) (also known as CD13) as a binding receptor allows infection of a diverse range of species 13,20 . Lednicky et al., identified PDCoV strains in plasma samples of three Haitian children with acute undifferentiated febrile illness 12 confirming the ability of PDCoV to cause disease in humans.
SADS-CoV, a swine enteric alphacoronavirus, is a recent spillover from bats to pigs 1 . SADS is a highly pathogenic enteric CoV first reported in a fatal diarrhea outbreak in Guangdong province, China, in January 2017, causing the deaths of 24,693 newborn piglets 21 . This disease is caused by a novel strain of Rhinolophus bat coronavirus HKU2 1,17 . Zhou et al., provided significant evidence that the causative agent of SADS-CoV is a novel HKU2-related coronavirus that has 98.48% identity in genome sequence to HKU2 17 .
Transcriptomic analysis is a powerful application of Next Generation Sequencing (NGS) technology that allows the identification of pivotal genes and/or signaling pathways as well as biomarker and drug discovery for various diseases and novel therapeutics 22  As a result of the increasing availability of transcriptomic data, it is possible to identify different transcriptomic datasets under similar disease and control conditions that can help to elucidate novel pathways and genes with remarkable accuracy 22 . Therefore, transcriptomic analyses using different datasets are becoming increasingly useful. Krishnamoorthy  This study compares the transcriptomic profiles of epithelial cells infected with distant CoV such as: human intestinal epithelial cells (HIECs) infected with PDCoV 21 , normal human bronchial epithelial (NHBE) cells infected with SARS-CoV-2 41 , and porcine intestinal epithelial cells (IPEC-J2) infected with SADS-CoV 22 . We hypothesized that similar and unique aspects in the immune-associated response to coronavirus infection in epithelial cell lines exist. Comparison of the host response to highly pathogenic coronaviruses versus potential emerging human pathogen PDCoV will provide key knowledge in understanding and developing therapeutic targets for these diseases. To explore and compare the transcriptome profiles of epithelial cell lines infected by PDCoV, SARS-CoV-2, and SADS-CoV infection, this study was able to identify common differentially expressed genes (DEGs) and signaling pathways among phylogenetically distant CoV. To our knowledge, this is a novel comparative transcriptomic analysis of epithelial cells lines infected by coronaviruses differing at the genus level.

Results
Alpha-, beta-, and deltacoronavirus infections result in more differentially upregulated genes across 10 common immune response associated pathways. This  We found that the DEGs are associated with 10 common immune response associated pathways. The apoptosis signaling-, interferon signaling-, interleukin signaling-, T-cell activation-, TGF-β signaling-, and Ras signaling-pathways were mostly upregulated upon viral infection. Within these pathways, more genes were affected in the inflammation/cytokine signaling pathway in all CoVs in comparison to the other 9 pathways (Fig. 1, Fig. S1-S2). In the same way, more genes were affected in HIEC cells infected with PDCoV in this pathway and a small number of genes were affected in the Interferon and Jak-Stat signaling pathway in a majority of the cell lines ( Fig. 1, Fig. S8).
DEGs of each dataset were submitted to a gene set enrichment analysis (GSEA) using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment to identify pathways that were common upon viral infection, and then, manually represented using CorelDRAW2019 (Fig. 2). A total of 23 pathways were common and enriched in human and porcine cell lines at 24 hpi (Fig. 2). These included Cytokine-Cytokine Receptor interaction pathway (Fig. S1), Viral protein interaction pathway with cytokine and cytokine receptor (  Table 2).

There is a positive correlation between PDCoV and SARS-CoV-2 infections.
To see if there is a relationship between PDCoV and SARS-CoV-2 infections, we first identified genes that are upregulated and/or downregulated in both datasets (Fig. 4a). We focused on the analysis of DEGs that are upregulated and downregulated in HIEC and/or NHBE datasets. As a result, we were able to identify 2309 DEGs that were upregulated in HIEC cells infected with PDCoV and down-regulated in SARS-CoV-2 infected NHBE cells (Fig. 4b).
Similarly, we were able to identify 2567 DEGs that were upregulated in NHBE cells infected with SARS-CoV-2 and downregulated in HIEC cells infected with PDCoV (Fig. 4c). This study also enabled identification of common DEGs that are upregulated in HIEC cells infected with PDCoV and downregulated in NHBE cells infected with SARS-CoV-2 at 24 hpi across 9 pathways (Fig. S17) Fig. S17). Likewise, we were able to identify common DEGs that are upregulated in NHBE cells infected with PDCoV and downregulated in HIEC cells infected with SARS-CoV-2 at 24 hpi across 10 pathways (Fig. S18).  (Table S1d, Fig. S18). Next, we performed the correlation between the two datasets using Pearson's correlation coefficient and we found that there is a positive correlation between PDCoV and SARS-CoV-2 infections (Fig. S19).

Discussion
Novel coronaviruses crossing between animals and humans are likely to continue to be the source of future pandemics 68 . To anticipate therapeutics for future coronavirus outbreaks, there is an urgent need to understand potentially conserved therapeutic targets. This study was able to explore and compare the transcriptome profiles of epithelial cell lines infected by PDCoV, SARS-CoV-2, and SADS-CoV infection at 24 hpi; and to show that distant coronaviruses such as Alpha-, Beta-, and Deltacoronaviruses share common immune associated pathways and genes that could be used as conserved targets for drug discovery. Our analyses appear to be the first report of comparative transcriptome analysis of coronaviruses at the genus level, incorporating available public datasets  www.nature.com/scientificreports/ from alpha-, beta-, and deltacoronavirus host responses. It is important to point out that the main limitation of this research was lack of public transcriptome datasets available related to cell lines infected with different coronavirus groups. Future studies need to make their transcriptomic data readily available in order to make more accurate and complete comparisons. Through a gene level approach, we identified common DEGs in the transcription profile of each epithelial cell line infected with PDCoV, SARS-CoV-2, or SADS-CoV. These DEGs associated with 10 common immune related pathways that were mostly upregulated upon infection (Table 1, Fig. 1). These results are in line with publications reporting more than 80% of immune related genes being up-regulated 21,28,32 . Likewise, PDCoV, SARS-CoV-2, and SADS-CoV shared 21 DEGs across 10 common immune response associated pathways ( Table 2, Fig. 1, Fig. S15-S16). These genes are described per group/function as follow: MAPKAPK2, also known as MK2; MAPK10 and MAP3K8 are in the first group that corresponds to genes that are part of MAPK signaling and are involved in inflammatory processes, innate and adaptive immunity 42,43,46,58 . In this group, MAPK10 has recently been shown to be downregulated in patients with Down syndrome and COVID-19 46 . Other genes that are involved in inflammatory processes include ITGA4 (also known as CD49d), RELA, IRAK4, SOCS3, and genes that belong to the NF-κB signal transduction pathway such as NFKB1and NFKBIE 44,47,48,51,54,59,60,65 . In this case, NFKB1 is associated with the upregulation of inflammatory responses in patients with COVID-19 infection 47,48 .
PIK3CA, TRAF2, RELA, STAT2 and TICAM1 are part of a second group of genes that are involved in the activation of type I IFN expression and exhibit antiviral defense properties 45,48,[52][53][54]61,66 . In this group, RELA is involved in upregulation of inflammatory responses in patients with COVID-19 infection 48,54 . CDKN1A and BIRC3 are associated with another group of genes that are involved in apoptosis in some viral infections 49,57 .
Finally, some unique genes play key functions. For instance, PAK1 upregulation is associated with several cancers, malaria, influenza, HIV, and COVID-19 55 ; and RELB enhances viral transcription of some viruses in the nucleus 67 (Table 2). Taken together, these results suggest that many of these genes are involved in inflammatory processes, activation of type I IFN expression, apoptosis; and several of them are involved in upregulation of inflammatory responses in patients with COVID-19 infection. We also speculate that, based on visualization of expression profiles of these 21 DEGs in a heatmap (Fig. 3), PDCoV and SADS-CoV are using similar strategies for efficient viral replication in the host cells despite species differences in comparison to SARS-CoV-2. This is expected since PDCoV and SADS-CoV are both swine-adapted enteric coronaviruses.
Consequently, a pathway level analysis identified pattern recognition receptor (PRR) signaling pathwayssuch as Toll-like receptors (TLR) (Fig. S4), NOD like receptors (NLR) (Fig. S5), and RIG-I-like receptors (RLRs) (Fig. S6)-upon coronavirus infections. Additionally, this study has found that PDCoV, SARS-CoV-2, and SADS-CoV shared 23 pathways that included Cytokine-Cytokine receptor interaction (Fig. S1), viral protein www.nature.com/scientificreports/ interaction with cytokine and cytokine receptors (Fig. S2), NF-Kappa B (Fig. S3), Cytosolic DNA signaling (Fig. S7), JAK-STAT (Fig. S8), IL-17 signaling (Fig. S9), TNF signaling (Fig. S10), Malaria signaling (Fig. S11), Influenza A (Fig. S12), Coronavirus Disease (Fig. S13), among others (Fig. 2). Some of these pathways, including the JAK-STAT signaling pathway and the cytosolic DNA sensing pathway, play a pivotal role in innate immune responses 69,70 . Pathways such as Cytokine-cytokine receptor interaction and Toll-like receptor signaling have top ranked functions correlated to gamma CoV infection in avian cells (Fig. 2) 34 . Moreover, pathways such as NF-Kappa B signaling, RIG-I like receptor signaling, NOD like receptor signaling, Cytosolic DNA sensing, Influenza A, Malaria, TNF-signaling, Cytokine-cytokine receptor interaction and Toll like receptor signaling were significantly enriched during the SADS-CoV infection 21 . GO terms associated with response to virus, response to cytokine, regulation of MAPK, programmed cell death and inflammatory response were enriched during SARS-CoV-2 infection 28 . Also, cytokine response GO terms, IL-17 signaling pathway, TNF signaling pathway, and NF-Kappa B signaling were the chief pathways associated with SARS-CoV-2 36 . In fact, NF-Kappa B has been shown to be involved in inflammatory responses upon infection with human coronaviruses. Similarly, acute respiratory failure in influenza and COVID-19 patients is associated with the upregulation of IL-6 (an interleukin that is part of Cytokine-Cytokine Receptor Interaction Pathway) 71,72 . In this study, IL-6 was also upregulated in HIEC cells infected with PDCoV and NHBE cells infected with SARS-CoV-2 (Fig. S1). Another study that compared SARS-CoV-2 with Ebola virus, H1N1 influenza virus, MERS-CoV and SARS-CoV found that genetic pathways associated with hepatitis B, the AGE-RAGE signaling system, malaria, influenza A and rheumatoid arthritis were the most significant pathways altered upon infection 36 . In our study, we also found enriched pathways related to malaria (Fig. 2, Fig. S11), hepatitis B (Fig. 2), and influenza A (Fig. 2, www.nature.com/scientificreports/ Fig. S12). These results indicate that pathways associated with innate immune responses are shared among different coronavirus genres upon viral infection. The differences between DEGs that are upregulated and downregulated in HIEC and/or NHBE datasets with PDCoV and SARS-CoV-2 ( Fig. 4) allow us to hypothesize that PDCoV must overcome different host-innate immune evasion strategies in human cells to be more successful in the progression of infection while SARS-CoV-2 has likely already overcome many of these barriers. We can also see this trend in the change of gene expression in different pathways that were enriched upon viral infection (Fig. S1-S13). These DEGs were mainly related to apoptosis signal and immune response regulation. Interestingly, furin was upregulated in NHBE cells infected with SARS-CoV-2 and down-regulated in HIEC cells infected with PDCoV. Furin has a noted role in viral pathogenesis; for coronaviruses, furin cleavage sites have been found widely in betacoronaviruses as well as in avian-origin gammacoronaviruses and certain canine and feline alphacoronaviruses (Table S1f.) 73 . MAPK3 was upregulated in NHBE cells infected with PDCoV and downregulated in HIEC cells infected with SARS-CoV-2 in almost all 10 pathways (Fig. S18 and Table S1d). Buggele et al. showed that MAPK3 and IRAK1 protein levels are reduced during influenza virus infection; while this gene, also known as ERK1 74 , can positively regulate RNA/protein synthesis in astroviruses, viral protein synthesis in alphaviruses, CVB3 replication, hepatitis C virus (HCV) genome synthesis, and avian leukosis virus replication and virus-induced tumorogenesis 75 . RAC1 was upregulated in NHBE cells infected with PDCoV and downregulated in HIEC cells infected with SARS-CoV-2 in B-cell, Inflammation, Ras and T-cell signaling pathway. RAC1 plays an important role in the regulation of virus entry, replication, and release [76][77][78][79] and the inhibition of this gene leads to enhanced virus production 76,80 . CALM3 that was also upregulated in NHBE cells infected with PDCoV and downregulated in HIEC cells infected with SARS-CoV-2, was found in B and T cell activation pathways. A study has found that inhibition of Calmodulin-Dependent Kinase Kinase blocks human cytomegalovirus and severely attenuates production of viral progeny 81 . Nevertheless, the response to PDCoV showed a positive correlation (Pearson correlation coefficient = 0.4625) to SARS-CoV-2 infection (Fig. S19). Positive correlation also indicates that DEGs in both datasets showed similar changes in expression and can also explain the large overlap of DEGs between these two viruses. Similar trends were found in the comparison between SARS-CoV-2 and respiratory syncytial virus (RSV) 82 . These results suggest that even though there is a positive correlation in the immune-associated response to PDCoV and SARS-CoV-2, these viruses could be using different strategies to improve viral fitness in the infected host. Further study and validation of this correlation is warranted.
In summary, host cells infected with members of the phylogenetically diverse Alpha-, Beta-, and Deltacoronavirus genres share common immune associated pathways and genes. Therapeutics modulating these pathways may be effective in treating current and future novel coronavirus outbreaks. Ten common immune associated pathways were mostly upregulated upon infection across phylogenetically divergent coronaviruses. These viruses shared 23 pathways and 21 DEGs. The 21 DEGs provide conserved targets for drug discovery and potential therapy against emerging coronaviruses. Finally, despite a positive correlation between PDCoV and SARS-CoV-2 infections, these viruses could be using different immune escape strategies. We speculate that PDCoV is currently less successful than SARS-CoV-2 at causing significant disease in humans due to a lack of human adaptation. However, there is always a possibility that PDCoV can evolve to adapt similar mechanisms that have been used by SARS-CoV-2. Incorporation of a deltacoronavirus into comparative transcriptome analysis provides an initial report of host responses during infection with coronaviruses at the genetically distant genus level.
Expression data pre-processing. Gene expression counts were identified from the alignment files in BAM format using Rsubread 85 . Raw count data was transformed to counts per million (CPM) and log-CPM using EdgeR 85 . Genes that were not expressed in any biologically significant levels (CPM less than 1) were discarded. Expression values were normalized using the trimmed mean of M-values (TMM) normalization 86 . Differential expression. Differential expression analysis was performed between infected and mock samples in each cell line. Briefly, dispersion of each gene was first estimated, followed by fitting generalized linear models (GLM) on the expression dataset. Within the R programming environment, a quasi-likelihood (QL) F-test method was used to test for differential expression 85 along with Limma and EdgeR statistical packages being utilized for analyses. Statistically significant DEGs between uninfected control and the infected cell lines were identified using a Bonferroni-Hochberg adjusted P value cut-off of 0.05 and log fold change of 1 (FDR < = 0.1). Correlation analysis. Correlation of expression of DEGs in mock and infected samples was performed for two datasets: HIEC cells infected with PDCoV and NHBE cells infected with SARS CoV-2. Briefly, correlation of absolute Log fold change of DEGs between the two datasets was performed with the use of the Pearson's correlation coefficient. Additionally, scatterplots depicting correlation of Log fold changes were generated. Both the correlation analysis and generation of scatterplots were performed using the R Statistical and Programming environment.

Data availability
The www.nature.com/scientificreports/