Transcriptomic analysis reveals crucial regulatory roles of immediate-early response genes and related signaling pathways in coronavirus infectious bronchitis virus infection.

Coronavirus infection of cells differentially regulates the expression of host genes and their related pathways. In this study, we present the transcriptomic profile of cells infected with gammacoronavirus infectious bronchitis virus (IBV). In IBV-infected human non-small cell lung carcinoma cells (H1299 cells), a total of 1162 differentially expressed genes (DEGs), including 984 upregulated and 178 downregulated genes, was identified. These DEGs were mainly enriched in MAPK and Wnt signaling pathways, and 5 out of the 10 top upregulated genes in all transcripts were immediate-early response genes (IEGs). In addition, the induction of 11 transcripts was validated in IBV-infected H1299 and Vero cells by RT-qPCR. The accuracy, reliability and genericity of the transcriptomic data were demonstrated by functional characterization of these IEGs in cells infected with different coronaviruses in our previous publications. This study provides a reliable transcriptomic profile of host genes and pathways regulated by coronavirus infection.

Introduction

Coronaviruses (CoVs), a family of enveloped, positive-sense and single-stranded RNA viruses, infect humans as well as many other mammals and birds, often causing intestinal, respiratory, neurological or systemic diseases differing in severity (Fung and Liu, 2019, 2021). CoVs are currently classified into four genera, i.e., alphacoronavirus (αCoV), betacoronavirus (βCoV), gammacoronavirus (γCoV) and deltacoronavirus (δCoV). Among them, αCoV and βCoV, such as SARS-CoV (Li et al., 2005), MERS-CoV (Ge et al., 2013) and SARS-COV-2 (Zhou et al., 2020), are originating from bats, while γCoV and δCoV from birds (Woo et al., 2012). CoVs can cross the species barrier and become the causative agents of deadly zoonotic diseases. Infectious bronchitis virus (IBV) is an avian γCoV and causes respiratory tract disease of chicken, seriously affecting the performance of both meat-type and egg-laying birds. IBV infection is one of the main factors leading to economic losses in the poultry industry (Liu et al., 2019). However, the pathogenesis of IBV and the host antiviral mechanisms against IBV infection have not been fully elucidated.

When a viral infection occurs, a series of host genes are induced in the infected cells to deal with stress responses, leading to dramatic changes at the transcriptional level of host genomes. Transcriptomic analysis of differential gene expression and differential splicing of mRNA by RNA-Seq technology is an indispensable tool for initial determination of such changes (Sahraeian et al., 2017; Stark and Grzelak and Hadfield, 2019). As a high-throughput sequencing method with the advantages of high sensitivity and repeatability (Sultan et al., 2008; Wang et al., 2009c), it enables to accurately measure gene expression in novel gene transfer transcripts, identify a number of variable or splicing genetic events, find or export this new gene transcription, and accurately analyze genes expressed in different genotypes and their different functions, revealing the dynamic changes of transcriptome in different time and space (Hitzemann et al., 2013; Wang et al., 2009c). A number of transcriptomic studies have been reported in IBV-infected chicken cells (Dinan et al., 2019; Lee et al., 2021), chicken embryos (Hashemi et al., 2021) and chicken organs (Cong et al., 2013; Hamzic et al., 2016; Hashemi et al., 2020; Yang et al., 2017) as well as African green monkey kidney-derived Vero cells (Liao et al., 2011), revealing numerous DEGs in various cellular pathways and processes during IBV infection. However, the functional involvement of these DEGs in IBV replication, pathogenesis and virus-host interactions are yet to be characterized, probably due to the lack of suitable regents and tools for chicken and monkey cells.

In this study, we report the transcriptomic analysis of DEGs in IBV-infected human H1299 cells. GO function and KEGG pathway analyses have shown that these DEGs and their pathways are mainly involved in cell proliferation, growth, differentiation, apoptosis and inflammatory response as well as a number of other cell functions. The functional involvement of several DEGs in the immediate-early response gene (IEG) families, as revealed by this transcriptomic analysis, has been studied in detail in a number of previous publications, unveiling crucial roles of these IEGs in regulation of the replication, pathogenesis and virus-host interaction of IBV and other coronaviruses.

Result

Illumina sequencing and quality control

To investigate the differential regulation of host gene expression in IBV-infected cells, H1299 cells were either infected with IBV-p65 strain or mock-treated with UV-inactivated IBV, harvested at 20 h post-infection for RNA extraction. Total RNAs prepared from three replicates of IBV-infected and mock-treated H1299 cells, respectively, were sequenced by Illumina HiSeq platform. 57.49 Gb Clean data were obtained in total after filtering, and the Clean data of each sample was more than 6.83 Gb. The Q30 base ratio of each sample was 92.60% or above (Table 1 ), suggesting that the Illumina sequencing data were of high quality and suitable for further analysis. Finally, genes with differential expression levels were identified, and the functional annotation and enrichment analysis were subsequently carried out.

Table 1

Statistics of sequencing data.

Samples	Clean reads	Clean bases	GC Content	%≥Q30
mockr1	38,601,277	11,560,332,282	51.19%	92.60%
mockr2	35,232,744	10,545,153,320	51.26%	93.12%
mockr3	35,124,755	10,515,265,664	51.85%	93.23%
IBVr1	22,828,237	6,836,051,330	49.81%	93.48%
IBVr2	31,304,033	9,374,728,182	50.01%	92.90%
IBVr3	30,437,471	9,117,547,640	49.62%	93.03%

Differentially expressed genes between mock-treated and IBV-infected cells

Transcriptome sequencing analysis was used to screen DEGs with fold change (FC) ≥ 2 and a false discovery rate (FDR) at an adjusted p-value < 0.05. As shown in Fig. 1 A, the distribution of up- and down-regulated genes is shown using a volcano plot and MA plot, based on comparison of gene expression levels between mock-treated and IBV-infected cells (Fig. 1A). The heat map of hierarchical clustering was also presented in Fig. 1B. A total of 1162 DEGs were identified, including 984 up-regulated and 178 down-regulated genes (Fig. 1C).

Fig. 1

Analysis of differentially expressed genes

A) Volcano plot and MA plot of differentially expressed genes (DEGs) between mock-treated and IBV-infected H1299 cells. Each dot in the Volcano and MA maps represents a gene. The green and red dots represent significantly down-regulated and up-regulated genes, respectively, while the black dots represent non-significant genes.

B) Heat map analysis of differentially expressed genes. Each group has three repetitions. Green or red colors indicate significantly down- and up-regulated genes, respectively.

C) Number of up-regulated and down-regulated genes between mock-treated and IBV-infected H1299 cells.

GO enrichment analysis of differentially expressed genes

All DEG transcripts were further functionally divided into three GO categories: biological process, molecular function and cellular component. Through GO function analysis, a total of 30 different GO terms were found, with ten GO terms in each of the three GO categories (Fig. 2 ). Among the biological process terms, the cellular process, biological regulation and single-organism process terms have the largest number of DEGs. In the molecular function category, binding, catalytic activity and nucleic acid transcription factor activity were the top three terms, and for the cellular component category, the top three GO terms were cell, cell part and organelle. The gene numbers in these three categories were 696, 696 and 552, respectively (Fig. 2).

Fig. 2

GO enrichment analysis of differentially expressed genes

The abscissa is the number of genes, the ordinate is GO classification.

Pathway enrichment analysis of differentially expressed genes

The functions of DEGs were then analyzed with KEGG classification. It showed that DEGs were assigned to cellular processes, environmental information processing, genetic information processing, human diseases and organismal systems (Fig. 3 ). Most DEGs are classified into the environmental information processing, human diseases and organismal systems. MAPK signaling pathway contains the most DEGs in the environmental information processing, meanwhile, PI3K-AKT signaling pathway, RAS signaling pathway, TNF signaling pathway as well as cytokine-cytokine receptor interaction have the same number of DEGs. In the human diseases class, pathway in cancer, HTLV-I infection and transcriptional misregulation in cancer were the three top pathways with the greatest number of DEGs. “Osteoclast differentiation”, “oxytocin signaling pathway”, “Aon guidance”, and “estrogen signaling pathway” were the dominant terms in the organismal systems class. These results suggest that IBV infection may mainly affect the cell growth, differentiation and apoptosis. The top 3 enriched groups among the KEGG categories were pathway in cancer, HTLV-I infection and MAPK signaling pathway (Fig. 3).

Fig. 3

Pathway enrichment analysis of differentially expressed genes

The ordinate is the name of KEGG metabolic pathway, and the abscissa is the number of genes annotated to this pathway and their proportion to the total number of genes annotated.

KEGG pathway enrichment analysis was also performed, in which the top 20 pathways with the minimum of Q value were presented (Fig. 4 ). The main pathways activated by IBV infection were MAPK signaling pathway, Wnt signaling pathway, regulating stem cell pluripotency, osteoclast differentiation, cancer transcription disorder, tumor necrosis factor (TNF) signaling pathway and NF-κB signaling pathway. These pathways are mainly involved in cell functions such as proliferation, growth, differentiation, apoptosis, death and inflammatory response (Fig. 4). Activation of multiple signaling pathways in IBV-infected cells suggests that IBV may exploit these pathways to enhance replication and reproduction by regulating cell metabolism and other activities.

Fig. 4

KEGG pathway enrichment analysis

Each circle in the figure represents a KEGG pathway. The ordinate represents the pathway name, and the abscissa represents Enrichment Factor, which represents the ratio of the proportion of genes annotated to a pathway in the differential genes to the proportion of genes annotated to this pathway in all genes.

Statistics of DEGs in the significantly upregulated pathways and transcripts

To further study the functional roles of DEGs in IBV-infected cells, DEGs enriched in representative signaling pathways in transcriptomic analysis were analyzed. The results showed that cFOS, DUSP8, DUSP1, cJUN and NR4A1 were the most significantly expressed genes in the MAPK signaling pathway induced by IBV infection (Fig. 5 ). The top 5 genes in the Wnt signaling pathway were NFATC2, cJUN, WNT9A, WNT7B and PRKCG; in the signaling pathways regulating the pluripotency of stem cells were INHBA, ID4, WNT9A, KLF4 and LIF; in the pathway of transcriptional misregulation in cancer were FCGR1A, NR4A3, ETV7, DDIT3 and BIRC3; in the Circadian rhythm were BHLHE40, PER2, ARNTL, PER1 and CRY2; in the osteoclast differentiation were FOSB, cFOS, NFATC2, FCGR1A and cJUN; and in the TNF signaling pathway were cFOS, CXCL1, cJUN, PTGS2 and TNFAIP3 (Fig. 5).

Fig. 5

Statistics of DEGs in the significantly upregulated pathways and transcripts

X-axis represents the log value of multiple differences of gene expression levels between two samples, which is used to measure the difference in gene expression levels, and Y-axis is the genes name.

In all pathways, the top 10 genes with the most significant induction were FOXJ1, FOSB, cFOS, EGR4, CRB2, CASS4, ZNF474, RRAD, EGR3 and EGR1. These genes play important roles in cell activity, differentiation, apoptosis, inflammation and other pathological processes. Taken together, IBV infection activates the MAPK signaling pathway, Wnt signaling pathway, transcriptional misregulation in cancer, osteoclast differentiation, TNF signaling pathway and other important cell signaling pathways, regulating cell growth, differentiation, apoptosis, death and other physiological processes.

Validation for the expression of differential genes by RT-qPCR and Western blot

To verify the transcriptome sequencing results, time-course experiments were carried out in IBV-infected H1299 and Vero cells. The infected cells were collected at 0, 4, 8, 12, 16, 20 and 24 h post-infection, respectively, and the mock-treated cells were collected at 24 h post-treatment for RT-qPCR verification. Eleven up-regulated genes, CPEB3, cJUN, JUNB, JUND, cFOS, FOSB, SGK1, GADD45A, GADD45B, RGS2 and RND3, were selected for the verification and GAPDH was used as an internal reference gene for qPCR calculation and processing. As shown in Fig. 6 , the expression levels of these 11 genes increased gradually over time in IBV-infected H1299 cells, with some genes decreased after reaching the peak value, while others still showed an upward trend at the final time point (24 h). Among them, cJUN, JUNB, JUND, cFOS and FOSB genes belong to the AP-1 family, and their up-regulation levels significantly increased, especially cFOS and FOSB genes (Fig. 6A).

Fig. 6

Validation for the expression of differentially regulated genes in H1299 and Vero cells

A) Cells were infected with IBV at an MOI of approximately 2, and harvested at indicated time points for RNA extraction. Equal amounts of total RNA were reverse-transcribed. The level of IBV genomic RNA (IBV gRNA) and mRNA levels of DEGs were determined by quantitative PCR.

B) H1299 cells were infected with IBV (MOI∼2), or mock-treated. Cell lysates were harvested at 16 h and 20 h post-infection, respectively, and subjected to Western blot analysis with indicated antibodies. Sizes of protein ladders in kDa were indicated on the left.

In IBV-infected Vero cells, the induction of these 11 genes was similar to that in IBV-infected H1299 cells (Fig. 6A). Their expression levels were also up-regulated with the increase of IBV replication, but were at relative lower levels compared with those in IBV-infected H1299 cells (Fig. 6A). Some genes, such as RGS2 and CPEB3, showed higher induction in IBV-infected Vero cells than did in IBV-infected H1299 cells, and a similar induction level of GADD45A was observed in both types of cells (Fig. 6A). Among these genes, cFOS and FOSB genes were also the highest induced in Vero cells (Fig. 6A). It was also noted that the expression of some of these 11 genes decreased after reaching the peak, while the expression of several genes still increased till the last sampling time point (Fig. 6A).

Furthermore, Western blot analysis of cFOS (representing FOS family genes cFOS and FOSB), cJUN (representing JUN family genes JUN, JUNB and JUND), SGK1, GADD45A (representing GADD45A and GADD45B), RGS2, CPEB3 and RND3 at the protein level was then conducted in IBV-infected H1299 cells. Significant induction of the expression of cFOS, JUN, SGK1 and GADD45A was detected in the infected H1299 cells at 16 and 20 hpi, respectively, compared with the mock control (Fig. 6B). However, only minor induction of the other three proteins was observed in the infected cells at the same time points (Fig. 6B). Taken together, these results validate the accuracy and reliability of transcriptome sequencing data.

Discussion

Virus infection of cells regulates the expression of many cellular genes and pathways, taking advantages of host cell machineries for the successful completion of viral life cycle. In this study, transcriptomic analysis of IBV-infected H1299 cells reveals that a variety of cell signaling pathways and a large number of host genes were differentially regulated by IBV infection. The top 10 genes with the most significant differential expression during IBV infection were FOXJ1, FOSB, cFOS, EGR4, CRB2, CASS4, ZNF474, RRAD, EGR3 and EGR1. Using this information as starting points, the regulatory roles of several DEGs and related cellular pathways in coronavirus replication and pathogenesis were studied in detail in our previous publications (Li et al., 2022; Yuan et al., 2020, 2022; Zhu et al., 2021).

In the top 10 upregulated genes, FOSB, cFOS, EGR4, EGR3 and EGR1 are IEGs, a class of genes that are rapidly activated after external stimulation (Healy et al., 2013; Lau and Nathans, 1987). These genes are mainly involved in normal cell growth and differentiation, as well as in intracellular information transmission and energy metabolism (Healy et al., 2013). The FOS family includes cFOS, FOSB, FOS-associated antigen 1 (FRA1) and FOS-associated antigen 2 (FRA2) (Tulchinsky, 2000). Although most FOS proteins exist in a single isoform, FOSB proteins exist in two isoforms, the full-length FOSB and shorter FOSB2 (Delta FOSB), both lack 101 amino acids at the carboxyl terminus (Dobrazanski et al., 1991; Nakabeppu and Nathans, 1991; Tulchinsky, 2000). All FOS proteins contain bZIP motifs, which include leucine zipper regions responsible for heterodimerization with Jun to form dimeric AP-1 complex (Tulchinsky, 2000). The AP-1 transcription factor plays a key role in cell cycle control, apoptosis, cell differentiation, oncogenic transformation and tumor progression (Tulchinsky, 2000). Transcriptomic analysis revealed AP-1 as a key regulator that participated in the tumor inhibition induced by green tea polyphenols (Pan et al., 2014). EGR4, EGR3 and EGR1, belonging to the early growth response (EGR) family and the zinc-finger transcription factor family, are induced by a wide variety of extracellular stimuli including activation, growth and differentiation signals, tissue injury, and apoptotic signals (Mookerjee-Basu et al., 2020). EGR transcription factors play a role in TGF-β-dependent profibrotic responses (Bhattacharyya et al., 2011; Kosla et al., 2013).

Other 5 genes in the top 10 upregulated genes are FOXJ1, CRB2, CASS4, ZNF474 and RRAD. FOXJ1, a transcription factor in the Forkhead Box family, plays a critical role in cilia formation of respiratory, reproductive and central nervous systems (Chen et al., 2013), is associated with autoimmunity and involved in the T cell activity and NF-kB pathway (Lin et al., 2004; Srivatsan and Peng, 2005). CRB2, the fission yeast checkpoint protein, is specifically involved in the DNA damage checkpoint essential for the activation of downstream effector kinase Chk1 (Du et al., 2003; Qu et al., 2012; Saka et al., 1997). CASS4, also named HEPL, is one of the proteins of crk-associated substrate (CAS) family that act as scaffolds to regulate protein complexes controlling migration and chemotaxis, apoptosis, cell cycle and differentiation (Tikhmyanova et al., 2010). ZNF474 is a zinc finger protein, but its molecular function is rarely reported. RRAD (Ras-related associated with diabetes), a small Ras-related GTPase, is frequently inactivated by DNA methylation of the CpG island in the promoter region causing its down-expression in cancer tissues (Wang et al., 2014). RRAD up-regulation may act as a negative feedback mechanism to counteract cellular senescence by reducing the level of reactive oxygen species (ROS) (Wei et al., 2019).

The main pathways activated by IBV infection include the MAPK and the Wnt signaling pathways. MAPK signaling pathway mediates the intracellular signal transduction related to a variety of cellular activities and functions, including cell activity, differentiation, survival, apoptosis, death and transformation (Dhillon et al., 2007; Flores et al., 2019; Keshet and Seger, 2010; Plotnikov et al., 2011; Sabio and Davis, 2014). MAPKs are divided into four groups, namely, ERK1/2, ERK5, p38 and JNK (Fung and Liu, 2019; Keshet and Seger, 2010), and can be activated in response to environmental stimuli (Fung and Liu, 2019; Keshet and Seger, 2010). In IBV-infected H1299 cells, cFOS, cJUN, nuclear receptor subfamily 4 group A (NR4A1), dual-specificity phosphatase 1 (DUSP1) and DUSP8 are the most significantly expressed genes in the MAPK signaling pathway. Among them, cFOS and cJUN can form a stable, heterodimeric complex AP-1, involving in cellular proliferation, transformation and death (Shaulian and Karin, 2002; Tulchinsky, 2000). DUSP8 and DUSP1 can specifically dephosphorylate the threonine and tyrosine residues on MAPKs, assisting in shuttling or anchoring MAPKs to control their activities (Jeffrey et al., 2007). DUSP8 was shown to control the basal and acute stress in adult cardiac myocytes and regulate cardiac ventricular remodeling by altering ERK1/2 signaling (Liu et al., 2016), whereas DUSP1 can inactivate ERK1/2, JNK1/2 and p38 (Auger-Messier et al., 2013; Guo et al., 2020). NR4A includes three members, NR4A1 (Nur77/NGFI-B), NR4A2 (Nurr1) and NR4A3 (NOR-1) (Safe et al., 2016). As an IEG, NR4A can be rapidly induced by a variety of intracellular and extracellular stimuli, regulating metabolism and immune response (Rodriguez-Calvo et al., 2017; Safe et al., 2016).

Wnt signaling pathway, including a canonical or Wnt/β-catenin-dependent pathway and a non-canonical or β-catenin-independent pathway, is an evolutionarily conserved signaling pathway that controls cell growth, differentiation, apoptosis and self-renewal (Anastas and Moon, 2013; Clevers, 2006). The top 5 genes in this pathway were NFATC2, cJUN, WNT9A, WNT7B and PRKCG. NFATC2, also known as NFAT1 or NFATp, is a member of the NFAT family; all members of this family are transcription factors containing a Rel homology domain, which can coordinate with AP-1 (JUN/FOS) to regulate gene transcription and coordinate effective immune response (Hogan et al., 2003).

The functional significance of several DEGs in regulating coronavirus replication and pathogenesis has been investigated based on this transcriptomic analysis. DUSP1 was shown to function as a negative regulator of the p38 MAPK to restrict the induction of IL-6 and IL-8 in IBV-infected cells (Liao et al., 2011), and more recently, to be involved in the MKK3-p38-MK2-ZFP36 axis in coronavirus infection-induced proinflammatory responses (Li et al., 2022). Systematic studies of the functional involvement of cFOS and cJUN in coronavirus infection revealed important roles of these AP-1 family genes in regulation of coronavirus infection-induced apoptosis and the expression of proinflammatory factors via the MAPK pathway (Fung and Liu, 2017; Yuan et al., 2020; Zhu et al., 2021). The growth arrest and DNA damage (GADD) family, including GADD34, GADD45α (GADD45A), GADD45β (GADD45B), GADD45γ and GADD153, were upregulated in IBV-infected cells, as validated by RT-qPCR. Among them, GADD34 and GADD153 were significantly upregulated by IBV infection-induced ER stress response and play important roles in regulating IBV replication and virus-host interaction (Liao et al., 2013; Wang et al., 2009b). Our studies also reveal the essential roles of EGR family genes, especially EGR1, in regulating coronavirus replication and pathogenesis, and SGK1 in virion assembly and release (Yuan et al., 2022; Zhu et al., in preparation).

In summary, transcriptomic analysis of IBV-infected H1299 cells identifies a total of 1162 DEGs, including 984 up-regulated and 178 down-regulated genes. These DEGs were shown to be enriched in various pathways, particularly in the MAPK signaling pathway. RT-qPCR validation confirms the induction of 11 upregulated genes in IBV-infected Vero and H1299 cells. More importantly, the accuracy, reliability and genericity of the transcriptomic data have been demonstrated in our previous studies, revealing the essential roles of a number of DEGs, especially several IEGs, in regulation of coronavirus replication, pathogenesis and virus-host interactions.

Materials and methods

Virus, cells and antibodies

The egg-adapted Beaudette strain of IBV (ATCC VR-22) was obtained from the American Type Culture Collection (ATCC) and adapted to Vero cells as previously described (Fang et al., 2005; Lim and Liu, 1998; Wang et al., 2009a). H1299 cell line is a human non-small cell lung carcinoma cell line derived from the lymph node, and Vero cells were isolated from kidney epithelial cells extracted from an African green monkey (Liao et al., 2011). H1299 and Vero cells were cultured in Dulbecco's modified Eagle's medium (DMEM, Gibco) supplemented with 6% fetal bovine serum (FBS), 100 U/ml penicillin, and 100 g/ml streptomycin. Antibodies against cFOS (#2250), cJUN (#9165), SGK1 (#12103), GADD45A (#4632), β-actin (#4967) were purchased from Cell Signaling Technology. RGS2 (10,678-1AP) and RND3 (66228-1-Ig) were purchased from Proteintech. CPEB3 (A15402) were purchased from ABclonal. Goat anti-rabbit IgG H&L (Alexa Fluor® 488) (ab150077) was purchased from Abcam. IBV N protein was used in this study as previously described (Fung and Liu, 2017; Yuan et al., 2020).

RNA extraction

Total RNA was extracted using the TRIzol reagent (Invitrogen) according to the manufacturer's instructions. Briefly, cells were lysed with 1 ml TRIzol per 10 cm2 effective growth area, and the lysates were vigorously mixed with one-fifth volume of chloroform. The mixture was then centrifuged at 12,000g at 4 °C for 15 min, and the aqueous phase was mixed with an equal volume of isopropanol. The RNA was precipitated by centrifugation at 12,000g at 4 °C for 15 min, washed twice with 70% ethanol, and dissolved in 30–50 μl RNase-free water.

Transcriptomic analysis

Sample collection and preparation

RNA quantification and qualification: RNA degradation and contamination were monitored on 1% agarose gels, RNA purity was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA), RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, USA), and RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation for Transcriptome sequencing: 1 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-), and second strand cDNA synthesis was 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, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 200–250 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA), 3 μl USER Enzyme (NEB, USA) was then used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer, and PCR products were purified (AMPure XP system). The library quality was assessed on the Agilent Bioanalyzer 2100 system.

Clustering and sequencing: clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v4-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and paired-end reads were generated.

Data analysis

Quality control: raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.

Comparative analysis: the adaptor sequences and low-quality sequence reads were removed from the data sets. Raw sequences were transformed into clean reads after data processing. These clean reads were then mapped to the reference genome sequence. Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome. Tophat2 tools soft were used to map with reference genome.

Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (A manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database) and GO (Gene Ontology).

Quantification of gene expression levels: quantification of gene expression levels was estimated by fragments per kilobase of transcript per million fragments mapped. The formula is shown as follows: FPKM=(cDNA Fragments)/[Mapped Fragments (millions) ⅹ Transcript Length (kb)].

Differential expression analysis: differential expression analysis of two conditions/groups was performed using the DESeq R package (1.10.1). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 found by DESeq were assigned as differentially expressed.

GO enrichment analysis: Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseq R packages based Wallenius non-central hyper-geometric distribution (Young et al., 2010), which can adjust for gene length bias in DEGs.

KEGG pathway enrichment analysis: KEGG (Kanehisa et al., 2008) is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS (Mao et al., 2005) software to test the statistical enrichment of differential expression genes in KEGG pathways.

RT-qPCR and Western blot analysis

Total RNA was reverse-transcribed using the FastKing gDNA Dispelling RT SuperMix kit (Tiangen) according to the manufacturer's instructions. Briefly, 1 μg total RNA was mixed with 2 μl 5 x FastKing-RT SuperMix (containing RT enzyme, RNase inhibitor, random primers, oligo dT primer, dNTP and reaction buffer) in a 10 μl reaction mixture. Using a thermocycler, reverse transcription was performed at 42 °C for 15 min and the RT enzyme was then inactivated at 95 °C for 3 min. The cDNA was then diluted 20-fold with RNase-free water for quantitative PCR (qPCR) analysis, using the Talent qPCR PreMix SYBR Green kit (Tiangen) according to the manufacturer's instructions. Briefly, 8.4 μl diluted cDNA was mixed with 10 μl 2x qPCR PreMix, 0.4 μl 50x ROX, 0.6 μl 10 μM forward primer, and 0.6 μl 10 μM reverse primer for a 20 μl reaction mixture. The qPCR analysis was performed using a QuantStudio 3 Real-Time PCR System (Applied Biosystems). The standard protocol included enzyme activation at 50 °C for 3 min, initial denaturation at 95 °C for 3 min, followed by 40 cycles of denaturing (95 °C, 5 s) and annealing/extension (60 °C, 30 s) with fluorescent acquisition at the end of each cycle. The results obtained were in the form of cycle threshold (C ) values. Using the ΔΔC method, the relative abundance of a transcript was calculated using GAPDH as an internal control and normalized to the 0 h post-infection sample in the control group.

cDNAs were then subjected to qPCR using appropriate primers. The qPCR primers for IBV used in this study was listed as follows. IBV gRNA, 5′-GTTCTCGCATAAGGTCGGCTA-3′ and 5′-GCTCACTAAACACCACCAGAAC-3’; primer pairs for human and Vero cells used in the PCR were as follows. GAPDH, 5′-CTGGGCTACACTGAGCACC-3′ and 5′-AAGTGGTCGTTGAGGGCAATG-3’; FOSB, 5′-GCTGCAAGATCCCCTACGAAG-3′ and 5′-ACGAAGAAGTGTACGAAGGGTT-3’; JUNB, 5′-ACAAACTCCTGAAACCGAGCC-3′ and 5′-CGAGCCCTGACCAGAAAAGTA-3’; JUND, 5′-TCATCATCCAGTCCAACGGG-3′ and 5′-TTCTGCTTGTGTAAATCCTCCAG-3’; cFOS, 5′-GGGGCAAGGTGGAACAGTTAT-3′ and 5′-CCGCTTGGAGTGTATCAGTCA-3’; cJUN, 5′-AACAGGTGGCACAGCTTAAAC-3′ and 5′-CAACTGCTGCGTTAGCATGAG-3’; SGK1, 5′-GCAGAAGAAGTGTTCTATGCAGT-3′ and 5′-CCGCTCCGACATAATATGCTT-3’; GADD45A, 5′-GAGAGCAGAAGACCGAAAGGA-3′ and 5′-CACAACACCACGTTATCGGG-3’; GADD45B, 5′-TACGAGTCGGCCAAGTTGATG-3′ and 5′-GGATGAGCGTGAAGTGGATTT-3’; RGS2, 5′-AAGATTGGAAGACCCGTTTGAG-3′ and 5′-GCAAGACCATATTTGCTGGCT-3’; RND3, 5′-TTACACGGCCAGTTTTGAAATCG-3′ and 5′-GGGCGGACATTGTCATAGTAAG-3; CPEB3, 5′-CGCTGTCAAAATGGGGAACG-3′ and 5′-GTCCAAACCTGCGAAAGCTG-3’.

Western blot analysis was conducted as previously described (Yuan et al., 2020). In brief, cell lysates were harvested at the indicated time, and subjected to Western blot analysis using indicated antibodies. Gray values were calculated using Image J software.

Statistical analysis

The one-way ANOVA method was used to analyze the significant difference between the indicated sample and the respective control sample. Significance levels were presented by the p-value (ns, non-significant; *, p < 0.05; **, p < 0.01; ***, p < 0.001).

CRediT authorship contribution statement

Li Xia Yuan: Methodology, organized the study, did all of the experimental work, Formal analysis. Bei Yang: did all of the experimental work. To Sing Fung: Methodology, organized the study, Formal analysis. Rui Ai Chen: Methodology, organized the study. and.

Declaration of competing interest

The authors declare no conflict of interest.