Transcriptome analysis of primary chicken cells infected with infectious bronchitis virus strain K047-12 isolated in Korea.

Infectious bronchitis virus (IBV), an avian coronavirus, is highly contagious. Chickens with IBV infection develop acute pathogenesis in multiple organs, including the respiratory and urogenital tracts. Frequent recombination in the spike (S) glycoprotein gene has made vaccine strategies ineffective. To understand IBV pathogenesis, we analyzed the genetic distance between Korean IBV isolates and other coronaviruses, including SARS-CoV-2. To obtain comprehensive information about early immune responses such as innate cytokine production and associated immune regulation during IBV infection, we infected primary chicken embryonic kidney cells and performed transcriptome analysis. We observed that the functional pathways of innate immunity are regulated and confirmed expression of genes that coordinate early immune responses. Understanding the immune profile of the host cell may assist in vaccine development.

Introduction

Avian infectious bronchitis virus (IBV) is an important coronavirus that causes infectious bronchitis (IB) in the upper respiratory tract and spreads to remote organs [1, 2]. IBV is classified as a gammacoronavirus, while SARS-CoV-2 is classified as a betacoronavirus. Swine coronaviruses such as porcine epidemic diarrhea virus (PEDV) and porcine deltacoronavirus (PdCV) are classified as alpha- and deltacoronaviruses, respectively [3].

The spike (S) glycoprotein mediates entry of the coronavirus into the host cell. It is a major inducer of host immune responses and is a determinant of host range and cell tropism [3-5]. The S glycoprotein is cleaved into S1 and S2 subunits. The S1 protein facilitates viral attachment to the host cell receptor, and the S2 protein mediates the fusion between the viral and host cell membranes [4, 6]. Genetic modifications and recombination within the S1 coding region result in the generation of new viral subtypes. Consequently, vaccination is ineffective because genetic variations result in poor cross-protection [1, 2, 7].

When IBV exhibits renal tropism, it replicates in the tubules and induces interstitial nephritis and tubular pathology [8-10]. We previously reported the genetic features of a Korean IBV isolate (K047-12) that causes renal atrophy [11], and the current study was designed to analyze the host immune response to this strain. Using phylogenetic analysis of S1 protein sequences, we also examined the genetic relationship of IBV to other coronaviruses. Using chicken embryonic kidney (CEK) cells, we investigated changes in various immune pathways during acute infection with K047-12 by transcriptome analysis. Our results provide information about the role of host immune responses in renal pathogenesis at the cellular level when IBV affects kidney tissue. Moreover, understanding the interaction between avian coronaviruses and their hosts can contribute to preparation for the potential emergence of cross-species viral variants.

Materials and methods

The IBV K047-12 strain was isolated and propagated as reported previously [11]. Briefly, the virus was isolated in 2012 from the cecal tonsil of a 6-week-old layer in the city of Pocheon, Gyeonggi Province, Korea. The virus was transmitted via the in ovo route to 9- to 11-day-old specific-pathogen-free (SPF) chicken embryos. Before infecting CEK cells, the virus was passaged three times in SPF embryonated eggs (Sungmin Farm, Korea).

Six representative coronaviruses and three IBV strains were chosen for this study. The nucleotide sequences of the S glycoproteins were downloaded from the NCBI GenBank database as follows: HCoV-NL63 (NC_005831.2), TGEV (NC_038861.1), PEDV (NC_003436.1), SARS-CoV-2 (NC_045512.2), Bat-SL-CoVZXC21 (MG772934.1), IBV Beau-R (AJ311317.1), IBV M41-CK (MK728875.1), IBV K047-12 (MK618759.1), and PDCV (MN058072.1). S gene sequences were aligned by the Clustal W method using MEGA X software (version 10.2.2). A phylogenetic tree was constructed by the neighbor-joining method with 1,000 bootstrap replicates, using MEGA X software [12].

Next, CEK cells were prepared for transcriptome analysis by adapting a previously reported method [13]. Kidney tissues were collected from 18-day-old SPF eggs. Primary CEK cells were obtained by homogenization and seeded into 12-well plates (1 × 106 cells/well) in 1 ml of DMEM medium with 10% FBS. Cells were divided into a K047-12-infected group and an uninfected control group and cultured for 4-5 days until a monolayer was formed (~90% confluence). Cells were infected by replacing half of the medium with prewarmed medium containing an IBV virus suspension, and an equal amount of fresh medium without virus was used for the control cells. Cells were harvested at 48 hours postinfection and stored in TRIzol Reagent (Ambion, Austin, TX, USA) for extraction of total RNA. We repeated the infection experiments six to seven times, and the extracted RNA from a representative preparation was sent to Macrogen (Seoul, Korea) for RNA sequencing. The RNA was purified, fragmented, reverse transcribed, and amplified by PCR for sequencing. The expression levels of transcripts were calculated as the fragments per kb of transcript per million mapped reads (FPKM). The filtered data were normalized by trimmed mean of m-value (TMM) normalization using the calcNormFactors method in the edgeR R library [14-17]. The criteria for significant differential expression of genes between the two groups were as follows: fold change ≥ 1.5, raw p-value < 0.05. The statistical significance of differences in expression levels was determined using exactTest within edgeR. Gene ontology (GO) and functional annotation analysis for the DE genes was performed using gProfiler [18]. The significant data were analyzed using KEGG pathway enrichment. The transcriptome data for IBV strains (M41-CK, Beau-R, K047-12) were analyzed by hierarchical clustering (Euclidean distance and average linkage). The analysis was performed using Cluster 3.0 software and Java TreeView 1.2.0. The differential expression data for the M41-CK and Beau-R strains were obtained from a recent study [1].

For PCR, RNA was first reverse transcribed into cDNA using AccuPower RT mix (Bioneer, Daejeon, Korea) according to the manufacturer’s protocol. Quantitative PCR was then carried out in a 20-μL reaction mixture using Maxima SYBR Green/ROX qPCR Master Mix (Thermo Fisher, Waltham, MA, USA). All reactions were carried out in duplicate. The reaction conditions were 50°C for 2 min, 95°C for 10 min, and 40 amplification cycles consisting of 95°C for 15 s and 55°C (β-actin, IL-1β, IL-10, IL-12a, TLR4, TLR7, CCR2, CCR7) or 50°C (IL-6R, IL-12RH, IL-1RAP) for 1 min. The melt curve stage was performed at 95°C for 15 s, 60°C for 1 min, and 95°C for 15 s. Threshold cycle (Ct) values were determined, and relative gene expression levels were calculated by the 2-ΔΔCT method. mRNA levels were normalized against β-actin. Ten genes associated with immune responses were selected for validation at the mRNA level. Primers were designed based on published chicken genome sequences in the NCBI database. The primer sequences were obtained from published papers or designed using a primer design tool (Primer 3 and BLAST, NCBI) and are listed in Table 1.

Table 1

Primers used for qRT-PCR

Target	Forward primer (5'-3')	Reverse primer (3'-5')	Reference
β-actin	TTGTCCACCGCAAATGCTTC	AAGCCATGCCAATCTCGTCT	This study
Il-1β	ACAGAGATGGCGTTCGTTCCCGA	TCAGCTCGACGCTGTCGATGT	[49]
Il-10	CTGTCACCGCTTCTTCACCT	ATCAGCAGGTACTCCTCGAT	[50]
Il-12α	CAGAGCTGGGGAACCTCAAG	CATCTCTGCAGTGAGGGCAC	This study
Tlr4	TCTTTCAAGGTGCCACA	GGATATGCTTGTTTCCAC	This study
Tlr7	TCTGGACTTCTCTAACAACA	AATCTCATTCTCATTCATCATC	[51]
Ccr2	ATGCCAACAACAACGTTTGA	TGTTGCCTATGAAGCCAAA	[52]
Ccr7	CGGCTGAAGACCATGACAGA	CAGCCTGAACGATGGCAAAG	This study
Il6r	CGCCTGCTGGTGGAAGA	TTCACCCGGCAGACGAATTT	This study
Il2rg	CTACGTGCGCAGCAAGATCA	GACACTGTGAGGACACTGCGG	This study
Il1rap	TCGTGGATGCCAAAGTCACA	AGTGCTTGAGGCCGAATCTT	This study

Results

Sequences from the GenBank database were used to make a multiple alignment, and a phylogenetic tree was constructed to determine the genetic relationships between IBVs and other coronaviruses. We observed that the S1 protein sequences of IBV strains, including the Korean IBV isolate (K047-12), showed more similarity to those of SARS-CoV-2 and bat coronaviruses than to those of alpha- and deltacoronaviruses (Fig. 1A). Because the K047-12 strain causes renal damage, we infected primary CEK cells with the virus and performed RNA sequencing to investigate changes in expression of cellular genes. Out of 17,026 total genes examined, we identified 1,084 genes whose expression showed a >1.5-fold change. Of those genes, 787 were upregulated and 297 were downregulated (Fig. 1B).

Fig. 1

Phylogenetic analysis of IBV 047-12 and identification of differentially expressed genes during infection. (A) Neighbor-joining phylogenetic tree based on S1 sequences. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test is shown next to each branch. The tree is drawn to scale, with branch lengths in the same units as those of the genetic distances used to infer the phylogenetic tree. The IBV strains are highlighted. (B) Volcano plot showing p-values for comparisons of log 2-fold change in expression level between control cells and infected cells. The blue color represents downregulated genes, and the yellow bar represents upregulated genes (left). The number of differentially expressed genes is shown in the histogram (right). The criteria for significant differential expression were as follows: fold change ≥ 1.5, raw p-value < 0.05.

The differentially expressed genes were classified using GO enrichment analysis, which suggested significant changes in biological processes and molecular functions. In the analysis of biological processes, K047-12 infection significantly altered the expression of about 150 genes involved in leukocyte activation, inflammatory responses, and cytokine-mediated signaling pathways. Analysis of molecular functions revealed that pathways for cytokine and immune receptor activity were highly regulated (Fig. 2A). In addition, we examined functional classification by KEGG pathway analysis. In accordance with the GO enrichment analysis, the top-ranked functions were closely correlated with major immune responses, including metabolic pathways, cytokine-cytokine receptor interactions, and the Toll-like receptor (TLR) signaling pathway (Fig. 2B). The functional classification demonstrates that K047-12 infection of CEK cells activates the immune system.

Fig. 2

Functional classification of differentially expressed genes. (A) Classification by GO enrichment analysis. The top 10 terms are shown within biological processes and molecular function. (B) Heat map showing the results of gene-set enrichment analysis using a modified Fisher’s exact test on each pathway map. The legend shows p-values, where a p-value < 0.05 indicates a significantly enriched pathway term.

To identify strain-specific effects of the K047-12 strain on the host immune system, we compared our gene expression data with those obtained previously with other IBV strains. We retrieved the raw data for gene expression in primary renal cells that were infected with M41-CK or Beau-R in a previous study [1]. We then clustered genes by their functions in the immune response, cytokine responses, and cell activation. As shown in Fig. 1A, the S1 nucleotide sequence of M41-CK is more similar to that of Beau-R than to that of the K047-12 strain. However, the pattern of genes activated by K047-12 infection was more similar to that of M41-CK than to that of Beau-R (Fig. 3A). This suggests that the immune response is affected by factors other than the S1 sequence. Based on their functional classification, we identified highly regulated genes associated with innate immunity (Table 2). To verify the transcriptome analysis, we designed primers for 10 genes and confirmed their expression at the mRNA level. Following K047-12 infection, the levels of IL-1β, IL-10, IL-1RAP, and IL-2RG for inflammatory signaling, TLR4 and TLR7 for RNA-virus-specific pattern-recognition receptors (PRRs), and CCR2 and CCR7 for lymphocyte migration increased significantly (Fig. 3B). In contrast, IL-12a and IL-6R expression did not change significantly but did show an increasing trend. This result validates the transcriptome analysis and shows that renal cells are a source of critical immune regulators during the early stage of IBV infection.

Fig. 3

Classification of differentially expressed genes by immunologic function. (A) Hierarchical clustering of differential expression profiles of genes expressed in cells infected with strains K047-12, M41-CK, and Beau-R. Gene profiles were classified based on GO enrichment. The diagram indicates the correlation of the gene expression levels from all samples compared against each other. Colors represent the relative expression levels of genes, where yellow and blue indicate higher and lower expression, respectively (white color for no change). (B) Expression of genes involved in innate immunity, demonstrated by qRT-PCR. Cytokine and receptor-related molecules (upper) and TLRs and chemokine receptors (lower) expressed in primary chicken kidney cells are shown. Combined data from 6-7 experiments are shown, and error bars represent SEM. *, p ≤ 0.05 from uninfected groups by nonparametric and paired t-test

Table 2

Fold change in expression of highly upregulated genes induced by IBV infection

Gene	Description	Fold change(Infected / Control)
Il2rg	Interleukin 2 receptor subunit gamma	1.6447
Tf	Transferrin (ovotransferrin)	1.6469
Tlr4	Toll-like receptor 4	1.6504
Il6r	Interleukin 6 receptor	1.7668
Tlr15	Toll-like receptor 15	1.7790
Tlr7	Toll-like receptor 7	1.8052
Lcp1	Lymphocyte cytosolic protein 1	1.8258
Lbp	Lipopolysaccharide binding protein	1.8319
Il1rap	Interleukin 1 receptor accessory protein	1.8490
Stab1	Stabilin 1, transcript variant X2	1.8967
Il8l2	Interleukin 8-like 2	1.9669
Tlr2b	Toll-like receptor 2 family member B	1.9793
Cybb	Cytochrome b-245 beta chain	1.9834
Mmr1l4	Macrophage mannose receptor 1-like 4	2.2063
Cd48	CD48 molecule	2.2107
Ccl26	C-C motif chemokine ligand 26	2.2521
CclI8	Chemokine	2.3019
Il1r2	Interleukin 1 receptor type 2	2.3726
CclI7	Chemokine	2.3810
Cd80	CD80 molecule	2.4126
Cd72ag	CD72 antigen	3.0760
Ccr2	C-C motif chemokine receptor 2	3.1495
Marco	Macrophage receptor with collagenous structure	3.3576
Selp	Selectin P, transcript variant X2	3.9415
Ccl4	Chemokine (C-C motif) ligand 4	4.0789
Il1b	Interleukin 1, beta	4.2180
Ccr7	C-C motif chemokine receptor 7	4.2953
Il8l1	Interleukin 8-like 1	4.5621
Avd	Avidin	5.3144
Il4i1	Interleukin 4 induced 1	9.7158

Discussion

Complete prevention of IB by vaccination has not been very successful because of the rise of variant virus strains and unsuccessful vaccination strategies using inefficient routes, low vaccine titers, partial inactivation, and interference between vaccine strains [19, 20]. The respiratory track is the primary organ affected by IBV infection, but genetic recombination in the S1 coding region expands the range of remote tissues that can be infected [10, 21, 22]. Earlier studies on IBV focused more on viral activity and genetic features than on host-pathogen interactions [23-25]. Recently, we isolated an IBV strain that showed kidney tropism and reported its genetic characteristics by comparing IBV isolates from different regions [11]. In the current study, we examined the genetic relationship of the K047-12 strain to coronaviruses from different host species. Interestingly, IBV showed a relatively close relationship to betacoronaviruses such as SARS-CoV-2 and bat coronaviruses. Considering the repeated emergence of new variants amid the coronavirus pandemic (COVID-19), we speculate that a zoonotic coronavirus could arise through mutation of IBV. We therefore investigated how host cells respond to IBV infection and regulate their innate immune response. Infection by kidney-adapted IBV is harmful to chickens and results in systemic pathology. Thus, the transcriptome analysis of kidney tissue infected by an IBV strain with renal tropism is helpful for understanding the pathogenesis of infectious bronchitis [8, 26, 27]. As the first line of defense, innate immunity is activated within two days after antigenic stimulation [28]. At this initial stage of IBV infection, cytokines and chemokines are produced and co-stimulatory pathways are activated [29, 30]. To gain more precise information about the regulation of genes involved in innate immunity during kidney infection, we infected CEK cells with IBV isolate K047-12 for 48 hours and measured changes in the expression of genes associated with immune functions. The expression of 1,084 genes was found to change significantly; 787 genes were upregulated and 297 were downregulated. Assuming that the upregulated genes are linked to initial immune responses, we determined their functional significance using the GO enrichment and the KEGG pathway programs.

Functions for innate immunity regulation, such as leukocyte activation, inflammatory response, and cytokine-mediated signaling pathways, were activated by K047-12 infection. A recent study using the M41-CK strain showed a comparable gene regulation pattern [1]. For example, the RSAD2, SNX10, STAT1, TNIP2, IL-6, IL-8L1, and DLL1 genes, which are associated with cytokine pathways, were upregulated by infection with kidney-adapted IBV. Despite the genetic variability of the S1 gene, the hosts of different IBV strains might use common pathways for protection against the virus. However, it is debatable whether certain genetic modifications contribute significantly to pathogenicity and the host immune response. Although the modification of a few amino acids is enough to affect infectivity [31, 32], additional studies are necessary to understand the details of immune regulation. Moreover, our examination suggests the lack of a clear relationship between variations in the S1 protein and the host response. In a phylogenetic analysis of S1 protein sequences, the genetic distances between K047-12 and the M41-CK strain and between K047-12 and the Beau-R strain were similar (Fig. 1A) However, the immune response to K047-12 was more similar to the response to M41-CK than to the response to Beau-R (Fig. 3A). This observation emphasizes the advantage of transcriptome analysis in host cells over the viral sequence comparisons. Identification of major cytokine pathways that are regulated during acute IBV infection opens paths for determining how adaptive immunity manipulates systemic inflammation. We pinpointed several genes involved in innate immunity. In previous studies, infection with strain M41-CK activated IL-6 and IL-1β production in lung tissue, while an IBV strain with kidney tropism upregulated immune modulators such as IL-6, IL-10RA, IL-17RA, CCL4, and CCL17 [8, 33, 34]. Among the top 20 highly regulated genes, IL-1β, IL-10, IL-12a, IL-1RAP, IL-2RG, IL-6R, TLR4, TLR7, CCR2, and CCR7 were involved in all of the studies. We confirmed the expression of those genes using real-time PCR and identified IL-1β and IL-10 as signature genes. Produced by mesangial cells, IL-10 is a critical cytokine for pathophysiology in renal disorders. Despite its suppressive effect on inflammation, IL-10 can be upregulated during nephritis when inflammatory lymphocytes must be controlled and induction of regulatory T cells is required [35, 36]. Along with IL-6, IL-1β is secreted from renal mesangial cells and mediates inflammation. These are key activators for Th17 cell generation that develops in numerous renal disorders [37, 38].

During acute infection, TLR activation turns on intracellular signals to stimulate proinflammatory cytokines and their receptors, followed by generation of effector T cells, which is critical for protection, but excessive cytokine secretion is detrimental to the host [39-41]. For example, TLR3 and TLR7 are critical mediators of RNA virus recognition and initiation of the innate immune response [33]. TLR3 and TLR4 are also activated in SARS coronavirus infection [42, 43]. Considering that TLR4 and TLR7 are upregulated by K047-12 infection, we speculate that IBV might share common features with betacoronaviruses in the regulation of PRR recognition. Existing knowledge from rodent and human immunology about the regulation of chemokine receptors was applied to chickens. CCR2 receives the CCL2 signal and activates monocyte migration for viral clearance, and its deficiency is therefore fatal during infection [44]. CCR7 is activated by CCL19 and CCL21, and its critical role in T cell migration has been demonstrated in chickens [45]. In addition, human patients with SARS coronavirus infections showed activation of CCR2 and CCR7 expression [46, 47]. We also observed CCR2 and CCR7 activation by K047-12, suggesting that the chicken immune system activates chemokine receptors analogous to those in mice and humans and directs immune cell migration and host protection.

Genetic modifications in the IBV genome can result in altered tissue tropism, which makes it difficult to save animals from endemic outbreaks. IBV strains with renal tropism are lethal because they attack both urogenital and systemic organs [48]. Therefore, to understand the immune responses in the target organs to viral infection, we performed a transcriptome analysis of primary renal cells and found that the genetic profile in immunologic functions for pathogen recognition, cytokine regulation, and signaling is closely associated with innate immune responses during viral infection. This study provides information about immune pathways that are activated by infection with avian viruses. Expanding our basic knowledge about avian immunology is important for vaccine development and characterization of new variant viruses.