Cell-Type-Specific Molecular Responses to Avian Reovirus Inoculation In Vitro

1. Introduction

The Avian reovirus (ARV) remains relevant as a pathogen of significant interest as it affects poultry worldwide [1], causing substantial economic losses to the poultry industry [2,3]. The virus can cause various clinical manifestations, such as viral arthritis and tenosynovitis [4,5], intestinal lesions [6], hepatitis [7,8], pancreatic lesions [6,9,10], myocarditis and hydropericardium [11,12,13], and immunosuppression [14,15,16,17].

While the ARV-induced pathogenesis, symptoms, and lesions have been extensively reported and discussed, the transcriptome-wide cellular response, particularly in primary chicken embryo cell cultures, has not been investigated. In primary chicken embryo fibroblasts (CEFs), a direct relationship between autophagy and ARV replication [18] and apoptosis following intra-endosomal virion disassembly [19] have been reported. In a secondary fibroblast cell line (DF-1), a sustained antiviral response mediated by interferon-stimulated genes (ISGs) was observed [20].

While the relevance of interferons (IFNs) in the context of ARV infection has been studied, their function remains ambiguous. ARV protein σA-mediated antagonism to IFNs in CEFs [21,22] and no significant upregulation of IFNs in DF-1 cells [20] have been indicated. Despite the viral protein-mediated resistance, induction of both type-α and type-β IFNs in CEFs [23], as well as in joints of specific-pathogen-free (SPF) chickens [24] has been indicated following infection with ARV. A recent study investigating the transcriptome of spleens of chickens infected with ARV suggested a possible downregulation of IFN-β by overexpression of interleukin-4-induced-1 [25]. Earlier investigations had demonstrated an in vivo induction of IFNs following intratracheal or subcutaneous inoculation of pathogenic or attenuated ARV strains into white leghorn chickens [26]. Later, the same authors reported a lack of induction of IFNs by four virulent ARV strains in chicken kidney (CK) and chicken embryo kidney (CEK) cells [27]. Contrastingly, ultraviolet-inactivated ARV induced IFNs effectively in these two cell types and attenuated ARV induced priming-dependent IFN response in aged CEFs [27]. These discrepancies in reports suggest that the ARV-induced transcriptional changes could be cell-type-dependent and warrant a more comprehensive analysis.

The objective of the current study was to investigate the transcriptome-wide changes in gene expression patterns at various timepoints post-inoculation in three different cell types, i.e. primary CEK and chicken embryo liver cells (CELi) and , as well as a macrophage-derived continuous cell line (HD11).

2. Materials and Methods

2.1. Cell Culture

Chicken embryo kidney and liver cells were prepared from specific pathogen-free eggs at 19 and 15 days of embryonation (DOE), respectively. Briefly, kidneys and livers were collected at respective DOE, minced with scissors, and homogenized in 0.05 % Trypsin-EDTA solution (Gibco, Grand Island, New York). The homogenate was filtered with a 40-micron filter, centrifuged at 100  g for 10 minutes, and dissolved in 50 mL of growth medium (composition described below). HD11 cells were kindly provided by Dr. Li Zhang (Mississippi State University, MS).

All three cell types were cultured in Dulbecco’s Modified Eagle Medium (DMEM) (Corning, Corning, NY) supplemented with 10% fetal bovine serum (FBS) (HyClone, Logan, UT) and 2% Penicillin-Streptomycin with L-Glutamine (Corning, Corning, NY) and 3% sodium pyruvate (Corning). The cells were maintained at 37°C in a humidified incubator with 5% CO₂.

2.2. S1133 Inoculum Preparation and Dose Determination

The ARV S1133 strain from a stock prepared previously [28] was passaged once in CELi cells. The harvested cell lysate was freeze-thawed three times to release the virus. The virus was titrated on CELi cells with 8 replicates per dilution, and the titer was calculated to be 10^9.33 TCID₅₀/mL, using the Reed and Muench method [29]. This supernatant was diluted in DMEM to prepare inoculum at the final dose of 10⁶ TCID₅₀/mL.

2.3. Virus Inoculation

The cells were grown in one six 6-well plates per cell type, each timepoint and treatment. The growth medium was removed, monolayers were washed with 1 mL of phosphate-buffered saline, and 100 μL of ARV S1133 at 10⁶ TCID₅₀/mL were inoculated onto the monolayers of each well. Control wells received the same volume of DMEM only. Immediately after inoculation of the virus, 1.5 mL of DMEM containing 2% FBS was added to each well for post-inoculation maintenance. At 6, 12, and 24 hours post-inoculation (hpi) each plate to be harvested was removed from the incubator one at a time, placed on ice, and the medium was removed. For immediate RNA release from the cells, 1 mL of RLT lysis buffer (Qiagen, Hilden, Germany) was added to each well. The cells were mixed with RLT using a 21-gauge needle syringe and pipetted into 1.5 mL tubes. The homogenate was frozen immediately at -80°C until RNA extraction.

2.4. Quantitative PCR for Viral Load

RNA was extracted from 5 of the 6 samples using the RNeasy Mini Kit (Qiagen, Hilden, Germany) as per the manufacturer’s instructions. Viral RNA load was quantified by targeting the M1 gene of ARV. Total RNA was denatured for 10 minutes at 95°C and RNA was reverse transcribed using the LunaScript RT SuperMix Kit (New England Biolabs, Ipswich, MA). The qPCR was performed using Forget-Me-Not™ Universal Probe qPCR Master Mix, with the following primers: forward (5’-ATG GCC TMT CTA GCC ACA CCT G-3’), reverse (5’-CAA CGA RAT RGC ATCA ATA GTAC-3’) and probe (5’-TGC TAG GAG TCG GTT CTC GTA-3’) [30]. GAPDH gene was targeted to normalize the viral loads, using the primers designed in-house as follows: forward (5’-TGG TGG CCA TCA ATG ATC CC-3’) and reverse (5’-ACC TGC ATC TGC CCA TTT GA-3’), and probe (5’-ACT GTC AAG GCT GAG AAC GG-3’). All PCRs were conducted using the qTOWER³ PCR Thermal Cycler (Analytik Jena, Jena, Germany). Amplification peaks were analyzed by qPCRSoft program version 4.1. Relative viral RNA loads were calculated using the formula: $\log \left(2^{-(\mathrm{G}\mathrm{A}\mathrm{P}\mathrm{D}\mathrm{H}\mathrm{ }{\mathrm{C}}_{\mathrm{t}}-\mathrm{A}\mathrm{R}\mathrm{V}\mathrm{ }{\mathrm{C}}_{\mathrm{t}}\mathrm{ })}\right)$.

2.5. RNA Sequencing

Total RNA was submitted to GENEWIZ (South Plainfield, NJ), where the samples were quantified using a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). RNA integrity was checked using Agilent TapeStation 4200 (Agilent Technologies, Palo Alto, CA, USA). ERCC RNA Spike-In Mix ThermoFisher Scientific, Waltham, MA) was added to normalized total RNA before library preparation following the manufacturer’s protocol.

The RNA sequencing libraries were prepared using the NEBNext Ultra II RNA Library Prep Kit for Illumina, following the manufacturer’s instructions (New England Biolabs, Ipswich, MA, USA). Briefly, mRNAs were initially enriched with Oligod(T) beads. Enriched mRNAs were fragmented by heating for 15 minutes at 94°C. First-strand and second-strand cDNA were subsequently synthesized. cDNA fragments were end-repaired and adenylated at 3’ends, and universal adapters were ligated to cDNA fragments, followed by ligation of oligonucleotides for indexing (barcoding) and library enrichment by PCR with limited cycles. The sequencing libraries were validated on the Agilent TapeStation (Agilent Technologies, Palo Alto, CA, USA), and quantified by using Qubit 2.0 Fluorometer (ThermoFisher Scientific, Waltham, MA, USA) as well as by quantitative PCR (KAPA Biosystems, Wilmington, MA, USA).

The sequencing libraries were clustered on one flow-cell lane. After clustering, the flow cell was loaded on the Illumina instrument (4000 or equivalent) according to the manufacturer’s instructions. The samples were sequenced using a 2x150-bp paired-end configuration. Image analysis and base calling were conducted by the Control software. Raw sequence data (.bcl files) generated from the sequencer were converted into fastq files and de-multiplexed using Illumina's bcl2fastq 2.17 software. One mismatch was allowed for index sequence identification.

2.6. Read Alignment, Principal Component and Differential Expression Analyses

Raw sequencing reads were processed using fastp program [31], version v0.23.4, for quality control and trimming. Reads were then aligned to the Gallus gallus genome (GenBank accession GCA_016699485.1) using HISAT2 [32] version 2.2.1. For quantification, featureCounts [33] version 2.0.1 was used to count the number of reads aligned against the reference. The principal component analysis (PCA) plots were generated using built-in stats package in R, and PERMANOVA test [34] was utilized to assess statistical significance. The differential expression analysis was performed on the output of featureCounts [33] subread module version 2.0.1 using edgeR [35] version 4.2.1. The counts of differentially expressed genes (DEGs) were filtered at log2 fold change > 1 and P < 0.05 adjusted with Benjamini-Hochberg (BH) correction [36] to account for the false discovery rate. The output was used to create bar plots for the numbers of DEGs identified for each cell type and time point. Genes detected in each cell type were pooled across timepoints and used to generate an UpSet plot [37] using the UpSetR library [38] version 1.4.0. The statistical analyses and plotting were performed using RStudio version 2024.04.2+764 [39,40]. The scripts used for all the analyses is available at (https://github.com/Zubair2021/ARV_Cell_Transcriptomics_2024).

2.7. Pathway Annotation and Comparison

The lists of DEGs obtained from edgeR output were used as input to compare pathways enriched in each cell type using Metascape [41], version 3.5.20240901. The pathway networks were visually enhanced using Cytoscape version 3.10.1 [42], and the most significant pathway of the cluster was annotated using AutoAnnotate [43] version 1.3.0.

2.8. Protein-protein Interaction (PPI) Analysis

The interactions of proteins for each cell type were analyzed by inputting the lists of DEGs into the STRING database [44]. A PPI clustering pattern was observed only for CELi cells and is reported herein. The cluster networks were visually enhanced using Cytoscape version 3.10.1 [42].

2.9. Isoform Switch Analysis

Quality trimmed reads were aligned using splice-aware mapper STAR [45] version 2.7.11b, followed by transcript assembly with StringTie [46] version 2.2.3. Ballgown-type outputs [47] were generated comprising transcript abundance estimates, and differential isoform expression was conducted with IsoformSwitchAnalyzeR [48] version 2.4.0 to predict isoform usage following ARV inoculation.

3. Results

3.1. Avian Reovirus Replication in Different Cell Types

Avian reovirus S1133 replication was evaluated across the three different cell types, CEK, CELi, and HD11, at multiple time points post-inoculation. Figure 1 shows the relative viral RNA levels, quantified using qPCR targeting the ARV M1 gene and normalized against the housekeeping gene GAPDH. The data revealed a time-dependent increase in viral RNA across all cell types, with the highest levels observed in CELi cells at 24 hpi. The replication was greater in CEK cells compared to HD11 cells.

3.2. Principal Component Analysis of Gene Expression

Principal component analysis was performed to explore the variance in gene expression profiles between infected and control samples across time points and cell types. The PCA plot (Figure 2) revealed distinct clustering based on cell type. Among CEK and HD11 samples, some clustering was observed based on timepoints (Figure 3).

The statistical analysis using PERMANOVA test revealed statistically significant (P < 0.001) Euclidean distance between cell types suggesting that the variance in gene expression was attributed to cell type differences. While CEK and HD11 samples had a time-dependent distribution of samples, the PERMANOVA showed these were not statistically significant (P = 0.307). The treatment had the least effect on the clustering patterns and did not achieve statistical significance (P = 0.401).

3.3. Differential Gene Expression Analysis

Differential gene expression compared to uninoculated control cells was assessed across the three cell types infected with ARV S1133 at 6, 12, and 24 hours post-inoculation. As summarized in the bar plot (Figure 4), CEK cells had 49 DEGs at 6 hpi, decreasing to 4 DEGs at 12 hpi, before rising again to 30 DEGs by 24 hpi. Similarly, CELi cells had 26, 13, and 31 DEGs at 6, 12, and 24 hpi, respectively.

Primary chicken embryo kidney (CEK) and liver (CELi) cells, and immortalized macrophage-like cells (HD11) were inoculated with either ARV S1133 or DMEM. In CEK and CELi cells, a decline in the number of DEGs was observed at 12 hpi. The highest number of DEGs was identified in HD11 cells, especially at 12 hpi.

In contrast, HD11 cells exhibited the highest DEG counts, with 57, 170, and 94 DEGs at 6, 12, and 24 hpi, respectively. As illustrated in Figure 5, a total of 82, 63, and 318 unique DEGs were identified for CEK, CELi, and HD11 respectively. Only one DEG, i.e., fatty acid-binding protein 4 (FABP4), was found to be common among all three cell types (Figure 5).

3.4. Protein-Protein Interaction Network

A PPI network was constructed to examine the interactions between genes involved in the cellular response to ARV S1133 at various timepoints post-inoculation. The analysis revealed key immune-related genes (Figure 6: interferon alpha-inducible protein 6 (IFI6) associated with inhibition of viral replication by promoting apoptosis; radical S-adenosyl methionine domain-containing 2 (RSAD2, viperin) involved in viral replication inhibition by disrupting lipid rafts necessary for virus budding; myxovirus resistance protein 1 (MX1), a GTPase linked to viral replication inhibition through interferons, and 2’-5’-oligoadenylate synthetase-like (OASL) which is involved in degradation of viral RNA by activating RNase, as key nodes in the network underscoring their importance in mediating the cellular defense against ARV S1133 infection

3.5. Pathway Enrichment: Gene Ontology Analysis

As illustrated in Figure 7, an analysis of unique DEGs for biological process enrichment indicated a distinct profile of significantly upregulated pathways for each cell type. In CEK cells, immune response-related pathways were significantly enriched, including lymphocyte chemotaxis (GO:0048247, P < 0.01) and cellular response to type-II interferon (GO:0071346, P < 0.01). Additional processes such as chemokine-mediated signaling and regulation of hydrolase activity were also enriched (P < 0.01).

In CELi cells, defense response pathways were strongly represented, with enrichment for defense response to symbiont (GO:0140546, P < 0.001) and external biotic stimulus (GO:0043207, P < 0.001). Response to virus and interspecies interaction processes were also significantly overrepresented (P < 0.001). In HD11 cells, coagulation-related processes were prominent, including blood coagulation (GO:0007596, P < 0.01) and negative regulation of endopeptidase activity (GO:0010951, P < 0.01). Similar processes like wound healing and regulation of body fluid levels were additionally enriched (P < 0.01).

3.6. Network Analysis of Pathways

The pathway analysis using DEGs from three cell types revealed distinct but interconnected pathways (Figure 8) that were significantly upregulated in different cell types following infection or treatment. The pie charts at each node indicate the variable percentage of DEGs from each cell type being representing enrichment of a certain pathway.

The clusters primarily associated with CELi cells displayed a significant enrichment in pathways related to viral replication, leukocyte differentiation, and host responses to viral invasion. Genes in this cluster are involved in processes such as viral genome replication, transcriptional regulation by viral proteins, and antiviral immune responses.

DEGs from CEK cells primarily contributed to pathways related to the complement and coagulation cascades. This indicates that these cells had enrichment of processes involved in coagulation and hemostasis, likely in response to inflammation or cell damage induced by infection. Pathways such as platelet degranulation and the response to elevated platelet cytosolic calcium levels were also prominent.

Interestingly, DEGs in HD11 cells contributed to pathways predominantly related to growth and metabolism. The pathways included axonogenesis and sensory organ development, among others.

3.7. Isoform Switching Analysis

Since CELi indicated upregulation of pathways related to defense against the virus, splicing isoform usage as investigated in the infected CELi cells. As shown in the volcano plot (Figure 9), significant isoform switches for transcripts from 16 genes (shown as red dots) were observed.

The top 10 exhibiting significant isoform switch based on q-value and their functions have been listed in Supplementary Table 1. Most of the genes showing changes in alternative RNA splicing patterns have been associated with cancer. One gene, ANKRD17, has been previously linked with antiviral immunity through RIG-1-like receptor-mediated signaling.

A notable switch in ANKRD17 isoform usage between control and ARV S1133-infected samples has been illustrated in Figure 10, where isoform XM_040694181.5 showed a substantial increase in expression following infection.

4. Discussion

The study aimed at exploring ARV-induced transcriptomic changes across three cell types CEK, CELi, and HD11 at 6, 12, and 24 hpi. The primary CEK and CELi cultures were selected because of their pathological relevance and heterogeneous cell populations [50,51], as they better represent gene expression patterns of a tissue compared to secondary cell lines with more homogeneous cell populations. The cell line HD11 was selected to ascertain the immune responses of macrophages due to their relevance in ARV dissemination across tissues [52], potential association with impaired phagocytosis [53,54] and immunosuppression [55,56].

The quantification of viral RNA levels at various timepoints demonstrated a cell-type-dependent replication efficiency of ARV S1133, with the highest viral loads observed in CELi cells, followed by CEK cells, and the least efficient replication in HD11 cells. A higher replication of ARV in CELi compared to CEK is concordant with the previous findings on differences in the sensitivity of these two cell types to ARV infection [57]. Moreover, these results corroborate the observations of lower ARV replication in kidneys compared to livers in vivo and hence, lesser induction of lesions in the former [58]. Conversely, a higher replication efficiency of infectious bronchitis virus (IBV) in CEK compared to CELi has been demonstrated previously [59], further affirming that replication efficiency of a virus in cultured cells could reflect its tissue tropism in vivo.

While ARV S1133 has been shown to replicate in peripheral blood monocyte cultures [60], and chicken bone marrow-derived macrophages as efficiently as in CK cells [61], a lower replication of S1133 in HD11 cells in our study showed a restricted efficiency of productive ARV infection in this secondary cell line [62]. Interestingly, despite exhibiting limited ARV replication, macrophage-like cell line had the highest number of DEGs at all time points. Such a negative correlation of viral loads with number of DEGs has been reported for swine macrophages infected with the African swine fever virus [63]. Contrastingly, this negative correlation of viral replication with host gene expression remains at odds with the observations on IBV-infected HD11 cells [64]. Since our approach included viral genomic RNA as well as mRNA quantification, the viral loads represent both viral genomic RNA copies as well as the viral mRNA, unlike minus-strand specific PCRs, which signify progeny virus genome replication [65].

Interestingly, a lower number of DEGs was observed at 12 hpi for CEK andCELi, with a subsequent increase again at 24 hpi. The balance between viral replication and host response in these cells could be a potential viral strategy to maintain a favorable environment for propagation while avoiding overwhelming host shut off [66,67]. Although neither CEK nor HD11 cells had any common DEGs across timepoints, three genes i.e. OASL, IFIT5, and IFI6, were consistently differentially regulated in CELi cells at all 3 timepoints tested. These interferon-stimulated genes (ISGs), along with others, were also part of the protein-protein interaction network generated using DEGs expressed in CELi cells at all timepoints. These genes appear to be a part of the core antiviral response against avian reovirus in vivo [24,25,68] and in vitro [20], as well as against other viruses [69,70,71,72,73,74]. Remarkably, despite the most efficient ARV replication in CELi cells, this cell type had the lowest counts of DEGs.

In the ARV-infected CEK cells, pathways pertaining to type-II IFNs and IL-1 were enriched, validating the findings of previous researchers [23,24,75]. This is in contrast with previous studies demonstrating ineffective induction of IFNs in CEK cells [27] and the downregulation of IFN-β in spleens of chickens infected with ARV [25]. The enrichment of pathways related to heterophil and lymphocyte chemotaxis and migration in CEK cells suggested chemokine-induced preferential recruitment of these immune cells [76].

Surprisingly, HD11 cells exhibited differential regulation of individual pathways associated with blood such as fibrinogenesis and complement cascade activation. Similar observations have been made with IBV-infected chicken kidneys, where DEGs related to complement factors and blood vessel-associated pathways were identified [77]. Moreover, the role of the complement system in immune responses has been described for the influenza virus [78], human immunodeficiency virus [79], Sindbis virus [80], and dengue virus [81] infections. Additionally, an enrichment in several pathways involved in the regulation of endopeptidase activity could be speculatively associated with reovirus replication. Cathepsins, a class of endopeptidases involved in lysosomal proteolysis, have been linked to reovirus entry and disassembly [82]. Since the outer capsid processing of various reoviruses in macrophage-like cells has been associated with endopeptidase cathepsin-S [83], the regulation of endopeptidase activity in ARV-infected HD11 cells suggests a potential cellular response to infection to prevent ARV entry and disassembly. These results contradict in vivo studies with other viruses where significantly increased endopeptidase activity was observed in the pancreas of reovirus type-3 infected suckling mice [84] and high endopeptidase expression was related to severe outcomes with SARS-CoV-2 infection [85].

The identification of fatty acid-binding protein 4 (FABP4) as the only common DEG among all three cell types is a novel finding in the context of ARV infection. While FABP4 has not been studied in the context of reovirus-induced arthritis or tenosynovitis, it has been described as a biomarker of human knee osteoarthritis, where patients exhibited significantly higher systemic and synovial FABP4 [86]. Interestingly, a pro-inflammatory role of FABP4 has been described in the pathogenesis of chronic tendinopathy in humans. The results were confirmed with experimental tendon degeneration in mice [87]. Moreover, a higher level of FABP4 in lungs and circulation of patients infected with SARS-CoV-2, and an experimental alleviation of viral loads, lung damage, and fibrosis in infected hamsters treated with FABP4 inhibitors underscored its importance as a biomarker [88]. Roles of FABP4 in cancer [89], neutrophil recruitment in Pseudomonas infections [90], and inflammatory gene expression in grass carp [91] have been described.

Effects on alternative splicing of host precursor mRNA, resulting in expression of different transcript isoforms and their ultimate translation into proteins with divergent functions, have been reported previously for a mammalian orthoreovirus [92], as well as for influenza viruses [93,94]. In the present study, a novel observation of significant isoform switches due to alternative splicing of host transcripts in infected CELi cells adds another layer of complexity to the current understanding of the host response against ARV. A notable switch in isoform of the ANKRD17 gene transcript was especially interesting, given the implication of the gene product in antiviral immune responses via the retinoic acid-inducible gene-I (RIG-I)-like receptor (RLR) signaling pathway in influenza virus infections [95,96], and pro-inflammatory responses in bacterial infections [97]. The human orthologues of some other DEGs exhibiting isoform switches in ARV-infected cells are associated with cellular response to cancer. Since the oncolytic potential of ARVs has been documented previously [98,99,100], the differential isoform usage among various cancer-linked genes found in the current study suggests that ARV-induced oncolysis could potentially be explained by further exploration of alternative splicing.

In conclusion, the study describes a detailed molecular snapshot of ARV infection in three cell types, where viral replication efficiency, host gene expression, and cellular responses remain cell-type-specific but somewhat interconnected and uniquely meaningful.