Transcriptional profiling of host gene expression in chicken embryo lung cells infected with laryngotracheitis virus

Gene expression profile of lung cells infected by ILTV

Primary chicken embryo lung cells at passage 1 were infected by the USDA reference strain of ILTV and cells were collected at 1, 3, 5, and 7 dpi. Cytopathic effects (CPE) were observed by 3 dpi, which became more severe by 5 dpi. Massive cell disruption was observed at 7 dpi (Figure 1). Total RNA was isolated from both controls and infected lung cells at each dpi time point and subjected to microarray analysis.

To control for dye bias effects, spike-in control mixtures were utilized by mixing with RNA samples according to the manufacturer's recommendations (See Methods). The spike-in RNA controls consisted of two sets of synthetic RNA mixtures derived from the Adenovirus E1A genes with different concentrations in each set [29]. The Agilent chicken 4 × 44 K oligo gene expression array contains 320 spike-in indicator spots to be hybridized with spike-in controls of both A mix, which was hybridized with Cy 3, and the B mix hybridized with Cy 5 on each array. These spike-in sets were mixed with either uninfected control or infected samples and co-hybridized to arrays. The ratio of signal intensities for all spike-in spots were calculated, evaluated, and revealed no significant dye effects on all array slides (data not shown) as reported previously [29]. All raw and normalized data were deposited in the Gene Expression Omnibus (GEO; accession number: GSE20630).

Normalized signal intensities were subjected to statistical analysis to find differentially expressed genes during ILTV infection in cultured embryonic lung cells. The 44 K array revealed 11,491 genes with significant signal intensities that were sorted by signal to noise ratio (SNR) >3, meaning that real (forward) signals of the samples were three times greater than background signals. In order to discover time course change in gene expression patterns, a model-based method [30] was used for clustering the gene expression profiles. A key drawback in heuristic clustering techniques is that it is difficult to determine the number of clusters a priori. The method enables the number of clusters to be determined by estimating the number of components in a multivariate normal mixture model from which the data are generated. The clustering analysis resulted in three gene groups (Figure 2). Group 1 included a total of 789 genes that showed significant differential expression in response to ILTV, Group 2 included 6,265 genes that displayed moderate alterations, and Group 3 included 4,437 genes that revealed no alterations during ILTV infection at four time points in chicken lung cells. Of the 789 genes in Group 1 exhibiting differential expression in response to ILTV (see Additional file 1), the top 10% (79 genes) were sorted by statistical analysis based on the highest value of standard deviations using the mean values of four different time points (Table 1). This approach highlights genes with more significant alterations in response to ILTV overtime. Out of the 789 genes, 390, 370, 320, and 422 genes were down-regulated, while 399, 419, 469, and 367 genes were up-regulated relative to uninfected cells at 1, 3, 5, and 7 dpi, respectively.

Quantitative reverse transcription-PCR (qRT-PCR)

Expression clustering

The pattern of differential gene expression over time can provide insights into biologically functional relevance among genes. In the present study, a model-based clustering method [30] was used to cluster alteration patterns for the 789 differentially expressed genes in response to ILTV infection and revealed 7 gene clusters exhibiting distinct expression patterns (Figure 3 and Additional file 2). The 287 genes placed in cluster (C) 1 showed only nominal increases at 3 and 5 dpi followed by decreased expression levels at 7 dpi that were similar to those at the onset of the experiment. The C2 representing 97 genes exhibited a dramatic increase in gene expression only at 7 dpi, whereas the expression levels of the 90 genes in C3 progressively declined at 5 and 7 dpi. Three genes in C4 showed higher expression during early infection (1 dpi), sharp increases at 3 and 5 dpi, followed by a slight decline at 7 dpi. Expression patterns of 9 genes in C5 showed slightly lower expression at 1 dpi relative to the other time points, a dramatically increase at 3 and 5 dpi, followed by decreased expression at 7 dpi. The 85 genes in C6 showed lower expression at 1 dpi followed by a progressive increase during the later time points, which was opposite to 218 genes in C7 that showed higher expression at 1 dpi followed by decreased expression at 3, 5, and 7 dpi. GenBank accession numbers for genes in each cluster are shown in the Additional file 2.

Interestingly, the genes in C4 that exhibited the highest expression during ILTV infection include cytokines (IL8 and IL1β) and a chemokine (CXC-K60), while in the C5, IL6 was most highly expressed. From these findings, it is reasonable to hypothesize that expression of functionally relevant genes such as cytokines and chemokines released in response to an immune challenge may be regulated similarly during a specific immunological challenge.

Functional Gene Ontology (GO)

Recently, new bioinformatics tools have been developed to facilitate efficient analysis of biological functionalities for large numbers of differentially expressed genes obtained from microarray analysis. By using the IPA program (http://www.ingenuity.com/), bioinformatics aspects of differentially expressed genes during ILTV infection were analyzed for the relevance of gene functionalities and gene networks. While 789 differentially expressed genes were used as the input number of genes, only 275 have been characterized with specific cellular functions according to the IPA program. Results obtained with the IPA program in terms of numbers of genes, biological functions of genes and categories were independently confirmed with a second bioinformatics tool, DAVID (D atabase for A nnotation, V isualization and I ntegrated D iscovery) version 6.7 (http://david.abcc.ncifcrf.gov/) (data not shown). The group of 275 differentially expressed genes was placed into 65 functional groups (see Additional file 3) and the top 25 functional groups of genes are displayed in Figure 4. The main categories for gene functionalities include diseases and disorders, molecular and cellular functions, and physiological system development. Gene information was repeatedly used in multiple groups of functionalities due to the multi-functional characteristics for designated genes. It can be seen that the cancer related function contains the highest number (140) of genes, while 125 genes were involved in genetic disorders, and 54 genes were grouped as inflammatory responses.

Gene network analysis

Gene network analysis, which represents the intermolecular connections among interacting genes based on functional knowledge inputs, was performed on the differentially expressed genes using IPA program (see Additional file 4). In this way, 21 possible gene networks for all days post ILTV infection were generated based on differential gene expression. Of these 21 gene networks, only 6 gene networks were identical at all time points following ILTV infection. Possibly, the reason why only 6 of 21 networks were identical might be because of differences in the sets of focus molecules which are generated from p-values and fold change values of differential gene expression that are used in IPA algorithms. The dynamics of alterations in gene expression for a subset of genes during the time course of ILTV infection can provide insights into biologically interacting genes within a network that display functional similarities. The most interactive network (network #1) is presented in Figure 5 whereas the other networks are shown in Additional file 5.

Network #1 is closely associated with a signaling pathway of IL6, which is a cytokine known to be involved in cell proliferation and inflammatory responses [32]. The top functions related with genes in network #1 involve cancer, gastrointestinal disease, and the cell cycle. Interestingly, expression of certain genes in network #1 such as CDC20 (cell division cycle 20 homolog), PTTG1 (pituitary tumor transforming 1), CDC2, and Cyclin B, which are associated with cellular proliferation in cell cycle progression, appeared to be inversely related to IL6 expression. The dynamics of alterations in gene expression over time during ILTV infection suggest that ILTV infection elevates IL6 expression followed by the inhibition of cellular proliferation. In contrast, expression patterns of HPGD (hydroxyprostaglandin dehydrogenase 15-NAD), SOCS (suppressors of cytokine signaling), JAK (Janus kinase 1), and NASP (nuclear autoantigenic sperm protein) were independent of the IL6 expression pattern (Figure 5A, B and Additional file 5A). JAK is known to enhance cellular proliferation through the signal transducer and activator of transcription (STAT) pathway that can be suppressed by IL6 signaling [33]. The consistent downregulation of JAK supports a role of JAK in the repression of cellular proliferation by ILTV infection. The top functions of genes in network #2 are involved with cellular compromise, connective tissue disorders, and post-translational modifications. Several heat shock proteins (HSP) were also focused in this network (see Additional file 5B). Heat shock proteins, especially the HSP70 family that serve as molecular chaperones, are known to interact with viral early immediate genes in HSV-1 genomic DNA replication [34]. Interestingly, since the expression of several HSPs in network #2 were downregulated, it is reasonable to hypothesize that the lower HSP through 5 dpi may lead to production of erroneous virion structures of ILTV that in turn results in low ILTV titers in tissue culture, which has been reported to barely exceed one infectious unit per cell [1, 35].

Network #3 contains genes for growth factors and matrix metalloproteinases (MMPs), and genes have top functions associated with endocrine system function and development, carbohydrate metabolism, and digestive system function and development (see Additional file 5C). Expression for growth factors and MMPs increased on 3 dpi and remain elevated through 7 dpi. This observation is in agreement with other reports that infection with herpesviruses, such as HSV and HCMV, lead to an increase in growth factor expression and MMPs for extracellular remodeling, tissue invasion and angiogenesis [36–38]. Networks #4 and #5 contain cytokine genes (IFNβ and IL1β), chemokine genes (CCL20 and CCL4), and genes of the NF-kB families (NF-kB and NFIB) with top functions that are involved in organism injury and abnormalities, antigen presentation, cell mediated immune responses, lipid metabolism, small molecule biochemistry, and molecular transport (see Additional file 5D and 5E). Finally, genes in network #6 contain IL1, NF-kB, and ID1 that function in cardiac inflammation, cardiovascular disease and in the inflammatory response (see Additional file 5F). Interacting molecules found in networks #4, 5, and 6 are mostly focused on the host immune responses against pathogenic infections.

The network analysis suggests that a large number of biological pathways, regulated by various sets of genes, closely interact with each other in host responsiveness during ILTV infection. More detailed interactions among genes showing altered expression levels in each network are currently under investigation to identify host-response mechanisms that may occur in conjunction with general immunological reactions during ILTV infections.

The fold changes in gene expression of key molecules associated with cellular immune response, cell signaling, MMP molecules, cytokines, chemokines, and cell proliferation were plotted individually (Figure 6). This was done since the molecules may help in clarifying the interaction of host lung cells with ILTV. Four matrix metalloproteinases (MMPs) including MMP 7, 13, 16, and 27 were differentially expressed during ILTV infections. For example, MMP 1, 2, and 9 were shown to function in cell invasions of primary human endothelial cell in Kaposi's sarcoma-associated herpesvirus (KSHV) pathogenesis [39]. Similarly the viral oncoprotein meq in MDV is known to activate MMP3 transcription [40]. Furthermore, the balance between MMP9 and tissue inhibitor of metalloproteinases 1 (TIMP1) was altered in human macrophages in HCMV infection, such that MMP9 activity declined in response to HCMV infection. Moreover, it was shown that HCMV infection may affect atherogenesis in mice through the control of MMP9 expression [38]. Taken together, these results suggest that MMPs generally play a role in herpesvirus pathogenesis, but different isoforms of MMPs may be capable of responding to specific herpesviruses. Furthermore, in the present study with chicken lung cells, expression of surfactant protein A (SFTPA-1; GenBank accession - AF411083; Table 1 and 3) was downregulated by ILTV infection. Since the SFTPA-1 is a transcriptional indicator of EGFR (epidermal growth hormone receptor) signaling pathway [41], the reduction in SFTPA-1 expression suggests that the EGFR signaling pathway is suppressed by ILTV infection. This conclusion is consistent with the reported downregulation of EGFR functions in HCMV infected human lung [41] and foreskin fibroblastic cells [42]. Additionally, our findings are similar to a report in which the mRNA expression of SFTPA-1 declined after inoculation of influenza A virus H9N2 into chicken lung cells [43].

Another molecule, ovotransferrin (Otrf; GenBank accession-X02009; Table 1 and 3) has both iron transport- and antibacterial-activities. In the present study, the expression level of Otrf increased in ILTV infected cells (Figure 6A). The finding is qualitatively similar to a possible role of Otrf in MDV infection that might protect the spread of MDV in chicken embryonic fibroblast cells [44]. Otrf also was found to accelerate the expression of immune response gene such as IFN-γ against MDV infection [45]. Furthermore, lactoferrin, which is a homologous form of transferrin in mammals, showed antiviral activity against canine herpesvirus [46]. It has been hypothesized that the anti-viral activity of Otrf may be similar to the anti-HSV capability of mammalian transferrins [47]. IL6 plays a role in both pro-inflammatory and anti-inflammatory responses [32]. The elevation of IL6 expression during virus infection has been well-characterized as one of the immune response activities in the pathogenesis of various viruses, such as Dengue virus in human [48], or HSV-1 in mice [49]. Also, the expression of IL6 is increased by KHSV-encoding miRNA [50]. Similarly, expression of IL8 and IL1β were upregulated in ILTV infected cells and also in MDV infected chicken lung [51]. These observations are consistent with previous findings of NF-kB activation via IL8 signaling pathway by human herpesvirus (HHV)-8 infection [52] and HHV-6 infection in peripheral blood mononuclear cell cultures [53], respectively.

In addition to the well-characterized host-virus interactions, a variety of unique responses detected with the microarray analysis of ILTV infection in the present study. For instance, the expression level of vasoactive intestinal peptide (VIP; GenBank accession number-U09350: Table 1 and 3) decreased dramatically in ILTV infection, but the functional consequences were not determined. Likewise, genes related to various metabolic enzyme functions such as acyl-CoA synthetase long chain-1 (ACSL1; GenBank Accession number-AJ851480; Table 1 and 3) were differentially expressed in ILTV infected cells, but the precise mechanisms in the host response have not been verified. Therefore, further investigations are being performed to identify unique and more deeply involved interactions between host lung cells and ILTV.

Cell culture and ILTV infection

Cell culture reagents were purchased from Invitrogen Life Technologies (Carlsbad, CA, USA). Chicken embryo lungs were isolated from 19 day old specific-pathogen free (SPF) chicken embryos (Charles River Laboratories, North Franklin, CT, USA). Lung tissues were homogenized and incubated in a 0.125% trypsin solution for 30 min at room temperature (25°C). Cells dissociated from lung tissues were suspended in a 1:1 ratio of mammary epithelial growth media (MEGM; Lonza, Rockland, ME, USA) and Dulbecco's Modified Eagle's Medium (DMEM, 0.45% glucose) plus 2% fetal bovine serum (FBS), 100 units/ml penicillin, 100 μg/ml streptomycin, and 2 mM L-glutamine in 10 cm tissue culture dishes (Sarstedt Inc., Newton, NC, USA) pretreated with 0.5% gelatin in PBS to improve cell adhesion. Cultured cells were grown at 39°C containing 5% CO₂ until cells reached confluent monolayers (2 to 4 days). The USDA reference strain of ILTV (National Veterinary Services Laboratories, Ames, IA, USA) was used to infect the chicken embryonic lung cells at a multiplicity of infection (m.o.i.) of 0.1. Infected cells were incubated at 37°C for 1 hr with rocking gently every 15 min. After the incubation, 10 ml of media, 1:1 MEGM/DMEM, were added to each culture dish, and the cells were incubated at 37°C in 5% CO₂ for up to 7 days. This research was performed under the permitted protocol approved by both the Institutional Biosafety Committee (IBC; permit number: 10007) of University of Arkansas and the Animal and Plant Health Inspection Service (APHIS; permit number: 102743) of United States Department of Agriculture (USDA).

Total RNA extraction

Total RNA was extracted from uninfected- or ILTV infected chicken embryonic lung cells at 1, 3, 5, and 7 dpi using TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) following the manufacturer's instructions. Total RNA was treated with DNase I (New England BioLabs Inc., Ipswich, MA, USA), and RNA was re-purified by the TRIzol reagent. The quality of RNA was checked by fractionation on an agarose gel (data not shown).

Probe labeling and microarray hybridization

A two color labeling microarray system was used to compare uninfected- and ILTV infected embryonic lung cells at 1, 3, 5, and 7 dpi. Fluorescently labeled complementary RNA (cRNA) probes were generated by using the Two Color Microarray Quick Labeling kit (Agilent Technologies, Palo Alto, CA, USA) and following the manufacturer's instructions. RNA spike-in controls were used to adjust possible dye effects following manufacturer's instructions. The Spike-in controls represent two sets of ten synthesized RNA mixtures derived from the Adenovirus E1A transcriptome with different concentrations in each set [29]. These spike-in sets were mixed with either uninfected control or infected samples and co-hybridized to arrays. Briefly, 2 μg of total RNA were mixed with Spike-ins and converted to cDNA using reverse transcriptase and oligo dT primers in which T7 promoter sequences were added. T7 RNA polymerase was used for the synthesis and labeling of cRNA with either Cy3 dye for the uninfected control or Cy5 dye for the ILTV infected samples. The fluorescently labeled cRNA probes were purified using the Qiagen RNeasy Mini Kit (Qiagen Inc., Valencia, CA, USA), and the concentration, fluorescent intensities, and quality of labeled cRNA probes were determined using a Nano-drop spectrophotometer (Thermo Scientific, Wilmington, DE, USA). An equal amount (825 ng) of Cy3 and Cy5 labeled cRNA probes were hybridized on a 4 × 44 K Agilent custom chicken oligo microarray (array ID: 017698). The hybridized slides were washed using a commercial kit package (Agilent Technologies, Palo Alto, CA, USA) and then scanned using a Genepix 4000B scanner (Molecular Devices Corporation, Sunnyvale, CA, USA) with the tolerance of saturation setting of 0.005%. Three biological replicates were conducted.

Microarray data collection and analysis

Background-corrected red and green intensities for each spot were used in the subsequent analysis. Global normalization based on local polynomial regression (loess) was applied to the intensities to remove effects that were due to undesirable systematic variations in microarray experiments rather than biological differences. The average values of the resulting normalized expression values in replicate hybridization sets were considered in the subsequent analysis. In order to identify a given set of genes that exhibited major alterations over time, a model-based clustering method [30] was employed, and the genes in the cluster were considered as differentially expressed over the time period. All analytic techniques were implemented in R (http://www.R-project.org).

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

Reverse transcription was performed with 3 μg of total RNA using Superscript II reverse transcriptase (Invitrogen Life Technologies, Carlsbad, CA, USA) with oligo dT_12-18 primers (Invitrogen Life Technologies, Carlsbad, CA, USA) following the manufacturer's instructions. The reverse-transcribed cDNA were diluted by 1:10 ratio and a portion (1 μl) of each product was subjected to qPCR under the following conditions: 40 cycles of 95°C for 30 s, gene-specific annealing temperature for 62°C for 1 min, extension for 30 s at 72°C, and a final extension at 72°C for 10 min. A non-template control and endogenous control (chicken GAPDH) were used for the relative quantification. The differential expression levels for the ILTV infected group were compared by the 2^-ΔΔCT method against the uninfected controls [31]. Primers for qRT-PCR were designed using Primer3 software (http://frodo.wi.mit.edu/cgi-bin/primer3/primer3.cgi) with these parameters: amplicon length, 95-100 bp; primer length, 18-27 nucleotides; primer melting temperature, 60-64°C; primer and amplicon GC content, 20-80%; difference in melting temperature between forward and reverse primers, 1-2°C. Primers were synthesized by Integrated DNA Technologies (Coralville, IA, USA). Primer information is listed in Table 2.

Bioinformatics

Functional interpretation of differentially expressed genes was analyzed in the context of gene ontology and molecular networks using the Ingenuity Pathways Analysis (IPA) 6.5 software (Ingenuity Systems^®; http://www.ingenuity.com). The differentially expressed genes were categorized, compared to genetic categories in the IPA database, and ranked according to p-values [54]. The IPA analysis determined the subcategories within each category which is supplied with an appropriate p-value and the number of genes identified. Since the size of the created network could potentially be enormous, the IPA software limited the number of molecules in the network to 35, leaving only the most important ones based on the number of connections for each focus gene (focus genes = a subset of uploaded significant genes having direct interactions with other genes in the database) to other significant genes [55].