Introduction

Human rhinovirus (HRV), a non-segmented positive sense RNA from Picornaviridae,, is the most prevalent human respiratory pathogen, consisted of three species, HRV-A, B and C. HRV is associated with common colds and concomitant sinusitis and otitis media episodes as well as lower respiratory tract infection (LRTI). Bronchiolitis and pneumonia among young children, exacerbations of asthma in older children, and exacerbations of chronic obstructive pulmonary disease (COPD) and pneumonia in older age are the predominant features associated with HRV at different age groups. Association of HRV with LRTIs has also been demonstrated in lower airway samples. Human airway epithelial cell is the primary target for HRV infection leading to interstitial pneumonia. An emerging body of evidence indicates that asthma exacerbations and COPD triggered by HRV alterations of pulmonary epithelial cells and subsequent airway inflammation rather than direct cytotoxicity on the cells.

Global transcriptional response of HRV-infected pulmonary epithelial cells has been investigated using in vitro and in vivo studies, however, the focus of the studies were mainly on a single species of HRV-A. Although clinico-epidemiological studies have revealed specific features associated with HRV species, the transcriptomic data is rather limited and have not been extensively understood with respect to different species. In the current study we hypothesized that HRV-B infection could deregulate expression of genes in a timely manner with potential to introduce novel targets attributed to species B for further investigations. To address this hypothesis, we performed gene expression analysis of the RNA samples obtained from alveolar epithelial cells infected with of HRV-B using Affymetrix Genechip technology. Our study reveals comprehensive gene expression profile of epithelial cell response to HRV-B infection across different time points in vitro. We have demonstrated remarkable gene expression changes in each time points post infection including novel virus-induced genes. Collectively, our data provides novel targets for further evaluation and to develop efficient anti-HRV therapeutic strategies.

Preparation of virus stock

Serotype 72 of HRV-B species was purchased from ATCC (ATCC VR-1182). H1HeLa cell line (ATCC CRL-1958) was grown in Leibovitz’s L-15 medium supplemented with 10% (v/v) fetal bovine serum (FBS), penicillin (IU/ml), streptomycin (100 μg/ml) and amphotericin B (2 μg/ml) and incubated at 37°C under CO2 free incubator. Stocks of HRV72 were propagated in H1Hela cells. H1Hela monolayer cells at 80% confluency were inoculated with HRV72 and virus absorption was carried out for 1.5 hours with gentle shaking. The cells were nourished with fresh media with 2% FBS and incubated at 34°C. The cell culture flasks were monitored daily and when HRV typical CPE covered more than 50% of the cells, they were harvested after three freeze-thaw cycles to rupture the membranes followed by centrifugation to pellet cell debris and stored at -80°C. Then the crude virus was concentrated using ultracentrifugation at 48k rpm for 2 hours using OptiSeal Polypropylene Tubes (optima L-100XP ultracentrifuge; Beckman Coulter, Inc.). The suspended impure virus was subjected to sucrose gradient using Beckman SW41 rotor to remove exogenous proteins and contaminants. The sucrose fractions were collected by peristaltic pump and analyzed using SDS-PAGE. SDS-PAGE analysis of sucrose fractions revealed that pure virus was present in high concentration at UP6 fraction as evident by the presence of three distinct bands corresponding to VP1-3 capsid proteins (S1 Fig). The fractions containing concentrated intact virus were dialyzed using 3500 MWCO Slide-A-Lyzer Dialysis cassettes G2 (3.5 ml) and subsequently followed by concentration step using concentrator 9K MWCO (Thermo Fisher Scientific Inc.) and aliquots were stored at -80°C.

Infectious center assay

H1Hela cells were plated over 6-well plates at a cell density of 4–5 × 105 cells per well and incubated overnight at 37°C. Tenfold virus suspensions were inoculated in duplicates and absorbed for 1.5 hrs. Then the inoculum was replaced with agarose overlay containing 1.5% of carboxymethylcellulose sodium salt (low viscosity; C5678 Sigma-Aldrich), magnesium chloride (0.04 M) and 2% FBS. When the layer was solidified the plates were turned upside-down and incubated for 5 days at 34°C. The cells were fixed with 3.7% formaldehyde and stained with 0.1% crystal violet solution. The number of plaques at each dilution was calculated as PFU. The titer of the purified virus was 5.2 × 109 PFU/mL, as depicted in S2 Fig.

Cell viability

A human adenocarcinomic alveolar basal epithelial cell line, A549, (ATCCCCL-185) was cultured in F12K medium (ATCC). A549 cells were seeded in 96-well plates at density of 12,500 cells/cm2 and incubated at 37°C in humidified incubator with 5% CO2 for overnight and then infected with various MOI (1, 5, 10, 15, 25, 30, 35 and 40) and incubated at 34°C. Cell viability was assessed at 24 and 48 hpi by adding 3-(4, 5-dimethylthiazol-2-yl)-2, 5-diphenyltetrazolium bromide (MTT) solution to each well and incubated for 4 hours in the dark condition. Then stopping solution was added to each well and incubated for 20 minutes. The absorbance was measured at 570 nm using a micro plate reader. Cell viability was also tested by Trypan blue exclusion method. A549 cells where seeded in 25 cm2 culture flasks and infected with virus MOI of 10 and 50. At 24 and 48 hpi the supernatant was centrifuged and the cells were trypsinized and collected from both infected and mock-infected flasks. The viable cells where counted under the phase contrast microscope using 0.4% (w/v) Trypan blue dye aqueous solution.

Infection of epithelial cells and total RNA preparation

The A549 cells grown in 25 cm2 culture flasks were infected with HRV72 at multiplicity of infection (MOI) of 10 or mock-infected at three independent replicates with a week interval. Total mRNA was extracted from both HRV72-infected and mock-infected A549 cells using GeneJET RNA purification kit (Thermo Fisher Scientific Inc., Germany), according to the manufacturer’s instructions. DNA contamination was removed using DNase I treatment (Worthington Biochemical Corporation). RNA concentration (at 260 nm) and purity (260/280 nm and 260/230 nm ratio) was measured by using a Thermo Scientific NanoDrop 1000 Spectrophotometer. The RNA integrity and DNA removal was confirmed by 1% denaturating agarose gel electrophoresis. The proportion of the full-length RNA was evaluated by microfluidic analysis using Agilent 2100 Bioanalysier using RNA LabChip Kit (Germany), as per manufacturer’s instructions. The samples with RNA Integrity Number (RIN) of 7 and greater were used for microarray hybridization.

Virus quantification

Viral RNA was extracted from cell culture supernatants at 6, 12, 24, and 48 hpi and also from third PBS wash of the cell monolayer before collecting the cells using GeneJET Viral DNA and RNA Purification Kit. Equal amounts of the purified RNA from both supernatant and cell lysate samples were reverse transcribed to obtain cDNA using Maxima First Strand cDNA Synthesis Kit (Thermo Fisher Scientific Inc., Germany). Universal Probe Library Assay design centre (https://lifescience.roche.com) was used to design primers to target a 94 nt region at IRES sequence (Table 1) to determine viral load in both samples using LightCycler480 Instrument II (Roche Life Science). Standard curve was generated by 5-fold dilutions of cDNA template of virus positive control and used to determine relative quantities of virus by extrapolating the Cp values of unknown samples against the linear regression. Expression of the target virus was normalized to the expression of the PSMB2 as housekeeping gene.

Microarray hybridization

Three independent replicates of RNA samples from mock-infected and HRV-infected A549-cells at MOI of 10 PFU were purified at 6, 12, 24, and 48 hpi and hybridized to microarray genechips for gene expression study. The samples were performed in Affymetrix GeneChipPrimeView Human Genome U133 Plus 2.0 microarray system includes more than 53,000 probes covering more than 20,000 genes and processed according to the manufacturer’s instructions (Affymetrix, Inc, USA). To achieve reproducible results, 100 ng of input RNA with RNA Integrity Number (RIN) of minimum 7.0 was applied as a standard amount for first-strand cDNA synthesis (at 42°C for 2 hours) primed with a T7-oligo (dT) primer using GeneChip 3’ IVT Express kit (P/N 901228, Affymetrix, USA). Second-strand cDNA was subsequently synthesized using DNA polymerase and RNase H (at 16°C for 1 hour). Linear RNA amplification using T7 transcription technology was employed in this kit. In the next step, the double-stranded cDNA was used as a template to amplify biotin-modified RNA (amplified RNA or aRNA). The aRNA was purified and concentration and purity was measured by using NanoDrop Spectrophotometer. To achieve optimal sensitivity of the assay, the labelled aRNA was fragmented to a mean size of 200 bases at 94°C and for 35 minutes before hybridization on gene chip array. Hybridization, staining and scanning of GeneChips were performed in Research Instruments Sdn. Bhd., Malaysia, using GenChip Hybridization, Wash, and Stain Kit (P/N 900720), according to the manufacturer’s instructions. Hybridization cocktail composed of 10 μg of fragmented aRNA and hybridization controls was hybridized to the GeneChip at 45°C for 16 hours and 60 rpm rotation in hybridization oven. Then the arrays were washed and stained in GeneChip Fluidics Station 450 prior to scanning with GeneChip Scanner 3000 (Affymetrix, Santa Clara, CA, USA).

Microarray data processing

The array fluorescent signals were obtained as DAT files and converted to probe intensity data (.CEL) files using Gene Chip Command Console (AGCC) Software (instrument control software available at: https://www.thermofisher.com/us/en/home/life-science/microarray-analysis/microarray-analysis-instruments-software-services/microarray-analysis-software.html). CEL files were normalized to remove effects related to systematic variations. Summarized expression values (CHP files) were created using Robust Multichip Analysis (RMA) summarization algorithm implemented in Affymetrix Expression Console Software 1.3 (available at: https://www.thermofisher.com/us/en/home/life-science/microarray-analysis/microarray-analysis-instruments-software-services/microarray-analysis-software.html). RMA workflow is based on quantile normalization of probe level signal intensities and has a general background correction to reduce interarray variability. The constructed CHP files were used for evaluation of individual hybridization success (QC control) and identification of outlier sample using Principle Component Analysis (PCA) in Affymetrix Expression Console Software 1.3. The quality of the starting RNA samples were evaluated using internal control housekeeping genes (GAPDH and ACTB). Entire target labelling process was monitored by adding four exogenous, premixed control spikes: lys, phe, thr, and dap. Efficiency of sample hybridization on the array was assessed by hybridization controls (bio B, bio C, bio D, and cre). QC analysis using Expression Console software showed that all the samples were located within threshold boundaries. The RNA samples hybridized to the array had RNA Integrity Number (RIN) of 7 or greater. The 3′/5′ ratios of internal control probe sets including GAPDH and ACTB were less than 3 for all the arrays indicating appropriate sample quality. The signal value profile was as expected with correct relative abundance for each of the spikes for both target labelling and sample hybridization indicating acceptable hybridization process. Normalized and log2 transformed data were used to identify differentially expressed genes (DEG) using one-way analysis of variance (ANOVA) and to construct heat-map dendrogram and hierarchical clustering in virus-infected, as compared to the mock-infected cells using Transcriptome Analysis Console v2.0 (available at: https://www.thermofisher.com/us/en/home/life-science/microarray-analysis/microarray-analysis-instruments-software-services/microarray-analysis-software.html). The data was filtered to select DEGs using a p-value of 0.05 and less with a ≥1.5-fold or ≤-1.5-fold cut-off. The microarray data have been submitted to Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) and assigned the accession number GSE87463.

Functional analysis

Functional Annotation Cluster tool provided by Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics resources 6.7, at NIAID/NIH, (available at: http://david.abcc.ncifcrf.gov/) was used for gene annotation and to identify functional categories enriched in differentially expressed genes at each time points. Modified Fisher’s exact test (EASE score of ≤ 0.01) was used to determine the probability of overexpressed gene ontology (GO) biological process terms. The tables were arranged according to enrichment score of ≥ 1.5.

Validation of microarray results by qRT-PCR

Representative differentially expressed genes from different functional groups including chemokine (CXCL8, CCL20, CXCL3), anti-apoptotic (BCL2A1), transcription factors (FOSL1, JUN, and EGR1) and signal transduction mediators (DUSP6 and GNB4) were selected to further validate by qRT-PCR. The expression profiles of the IFN genes including IFN-β, Il29 (IFN-λ1), and IL28 (IFN-λ2/λ3) were further evaluated using qRT-PCR. Universal Probe Library (Roche: https://lifescience.roche.com/) was used to design primers and probes (Table 1). In order to estimate PCR efficiency, standard curves were derived from 5×serial dilutions for all genes. The slop of standard curves was used to calculate PCR efficiency. The quantification was performed using LightCycler 480 Instrument II and LightCycler LC480 software version LCS480 1.5.0.39 under following settings: 95°C for 10 min denaturation followed by 45 cycles of 95°C for 10 sec, annealing at 60°C for 30 sec, 72°C for 30 sec. Comparative results were calculated and the expression of genes were normalized to the GAPDH in all samples at 6, 12, 24, and 48 hpi. The amplified products were further analyzed by agarose gel and also subjected to the melting curve analysis (65 to 95°C) at the end of the experiment to determine the specificity of the primers.

Statistical analysis

Data were analyzed using the Statistical Package for the Social Sciences (SPSS) version 18.0. All P-values were two-tailed and values of < 0.05 were considered as statistically significant. One-way analysis of variance (ANOVA) was applied to identify differentially expressed genes (DE genes). ANOVA with post hoc comparison using Tukey test was used to compare mean viral load between time points. Student’s independent sample t-test was applied to each data set for validation study by qRT-PCR. The data were expressed as the means ± SEM.

High viability of airway epithelial cells was noticed across the HRV infection

Clear and countable virus plaques were optimized in the presence of magnesium chloride. To comprehend whether HRV72 infection induces any cytotoxic effects, the viability of A549 cells was determined using MTT assay. The viability decreased as MOI increased and at MOI of 10 the cell viability was 95% at 6 and 12 hpi, and subsequently decreased to 86% and 69% at 24 and 48 hpi, respectively (S3 Fig). Trypan blue dye exclusion test also demonstrated high viability of the virus infected cells with 97% and 91% viability at 24 and 48 hpi, respectively. The A549 cell monolayers were inoculated with HRV72 at MOI of 10 PFU per cell, which was suggested as a peak viral titer in nasal samples of patients with HRV cold in another study. The time-course analysis of HRV72 infection in A549 cells was performed by applying this infection protocol. The virus induced specific cytopathic effect (CPE) characterized by shrinkage and rounding off the cells as well as detachment and destruction of cell layers, as observed by light microscopy (Fig 1). CPE started to be appeared at 12 hpi on A549 cells and continued to increase up to 48 hpi.

HRV increased in a time-dependent manner across the infection course

In order to confirm the active replication of the virus in A549 cells, the level of viral RNA was determined in samples collected from both cell lysate and cell culture supernatant at indicated time points. The cell lysate of the 1.5 hour post-infected monolayer cells was used as the baseline in virus quantification analysis. The third PBS wash was used as a baseline for virus quantification in supernatant. As shown in Fig 2. A, a significant increase of viral RNA in HRV72-infected A549 cell lysate was noticed at 12, 24, and 48 hpi compared to baseline at 1.5 hours of inoculation (p ≤ 0.001). Multiple comparisons also indicated significant differences among all time points (6, 12, 24, and 48 hpi) in any combination (p ≤ 0.001). The substantial fold change was started at 6 hpi (51.5 fold); continued at 48 hpi (190.4 fold) and highest viral load was observed at 12 hpi (457.6 fold). Analysis of the virus propagation in cell culture supernatant also demonstrated time-dependent increase of viral RNA up to 48 hpi (Fig 2B). Significant induction of viral RNA was observed at 24 and 48 hpi (p ≤ 0.001) with highest RNA fold at 48 hpi (295 Fold). According to the data acquired from experiments analyzing both intracellular and extracellular viral RNA level, it is apparent that HRV72 achieved higher viral loads during early time points after infection (highest at 12 hpi). This is in line with shorter incubation period in HRV infections compared to the other respiratory viruses.

HRV infection modified large number of genes at timely manner

The affymetrix microarray gene expression analysis was performed at four time points post infection because the time course analysis of the infection demonstrated that replication of this virus is time-dependent (Fig 2). A comparative analysis of three independent experiments between HRV-infected and mock-infected cells using a threshold of 1.5-fold identified collectively 991 differentially expressed (DE) genes across all time points, from which 459 (46%) were up-regulated and 532 (54%) were down-regulated. During infection down-regulated genes were surpassed up-regulated genes and the category of up-regulated genes consisted of largest fold changes. The samples were separated in two groups including virus-infected and mock-infected based on hierarchical clustering of gene expression patterns, as shown in Fig 3. A. Surprisingly, many genes were significantly affected at 6 hpi and continued to a less extent at 12 and 24 hpi but again increased substantially at 48 hpi (Fig 3B). Complete list of DE genes are provided in S1–S4 Tables. Using a more stringent 2-fold criterion the up- and down-regulated DE genes across four time pointes were combined and presented as separate gene list in Tables 2 and 3, respectively.

The overlaps of DE genes have been revealed across all-time points using Venn diagrams as shown in Fig 3C. During infection majority of the genes uniquely expressed in each time points (759 genes out of 991, 76.6%) and collectively less overlap of genes (23.4%) was found, which has been shown individually between time points as in Fig 3D (the complete list of overlap genes is presented in S5 Table). The minority of the affected genes at 6 hpi were common to other time points (50 out of 339 DE genes, 14.7%), with 41% (21 out of 51) in common with 12 hpi and 8% (4 out of 50) and 24% (12 out of 50) in overlap with 24 and 48 hpi, respectively and the rest were common in other combinations. Of 110 DE genes, 59 genes (53.6%) were overlapped including 35.6% (21 out of 59 DE genes) overlapped with 6hpi, and 7% (4 genes) and 17% (10) were common in 24 and 48 hpi, respectively. Of the 119 DE Genes, 64 genes (54%) uniquely expressed at this time points and the overlap genes are presented in Fig 3C and 3D. Taken together, the data revealed that HRV infection has significant effect on global gene expression and majority of DE genes were unique at each time point.

Diverse functional groups were enriched by HRV during the infection course

The DE genes with known and available functional annotations where focused for functional analysis using analysis Gene Ontology classification (DAVID). Diverse functional categories were enriched across the infection course are presented at Table 4 (p-value ≤ 0.5 and an enrichment score ≥ 1.5). Analysis of the DE genes at 6 hpi revealed significant enrichment of functional clusters that were associated with inflammatory response and chemokine activity. HRV72 infection resulted in the sudden induction of many chemokines and proinflammatory genes both in numbers and magnitude. HRV72 induced a mixture of CC (CCL2, CCL20) and CXC (CXCL1, CXCL3, CXCL5, CXCL8) chemokines. All of them showed maximum fold change at 6 hpi. Expression of CXCL8, CXCL3, and CXCL5 continued up to 12hpi and the rest of chemokines exclusively expressed at 6hpi (Fig 4A).

In contrast to the first time point which was limited to chemokine response, DAVID analysis at 12 hpi showed that DE genes fell into diverse functional categories such as regulation of apoptosis, anti-apoptosis, wound healing and blood vessel morphogenesis (Table 4). Combination of anti-apoptotic genes such as BCL2-related protein A1 (BCL2A1) and myeloid cell leukemia sequence 1, BCL2-related (MCL1) as well as pro-apoptotic genes such as tumor necrosis factor receptor superfamily, member 21(TNFRSF21) and member 19 (TNFRSF19), and tumor necrosis factor, alpha-induced protein 3 (TNFAIP3) were induced concurrently at 6 and 12 hpi (Fig 4B). Functional analysis revealed that limited groups were enriched at 24 and 48 hpi (Table 4). Genes related to the extracellular matrix (ECM) and remodeling were affected remarkably at 24 and 48 hpi. Expression of three fibrinogen chains, α, β and γ was down regulated at 6 hpi (-1.51 β, -1.52 α), 12 hpi (-2.08 α, -2.18 β, -2.09 γ), 24 hpi (-2.22 α, -2.64 β, and -2.29 γ) and 48 hpi (-2.01 α, -2.68 β, -1.85 γ). In spite of single induction of PLAU at 24 hpi, the result demonstrated a significant increase of PLAUR at 6, 12 and 24 hpi (1.91, 1.69, 1.62 fold, respectively). The expression of other ECM related genes were induced including thrombomoduline (THBD), matrix metallopeptidase (MMP).

It is of great importance to note that in first three time points of HRV infection, the expression of dual specificity phosphatase 6 (DUSP6; type II DUSP subgroup) was up-regulated by 2.86, 2.12 and 1.50 fold at 6, 12, and 24 hpi, respectively while the expression of DUSP22 (type I DUSP subgroup) was down-regulated at 6 hpi (-1.57 fold). Analysis showed that other signal transduction related genes were found to be significantly affected including GNB4, SOCS3, SOCS2, and PREX1 (Fig 4D). Early growth response 1 (EGR1) was the most strongly up-regulated gene across all the time points (25, 5, 1.5 fold at 12, 24 and 48 hpi, respectively). In addition to EGR1, a cluster of other transcription factors including FOS-like antigen 1 (FOSL1; 3.13 fold at 6hpi; 1.88 fold at 12hpi; 1.57 fold at 24hpi), FOS-like antigen 2 (FOSL2; 1.51 fold at 6hpi), FBJ murine osteosarcoma viral oncogene homolog (FOS; 2.94 fold at 12hpi), jun proto-oncogene (JUN; 3.45 at 12hpi; 1.66 at 24hpi), v-maf musculoaponeurotic fibrosarcoma oncogene homolog F (avian) (MAFF; 2.89 at 6hpi; 1.96 at 24hpi) and ETS-domain protein (SRF accessory protein 2) (ELK3; 1.97 fold at 6 and 1.52 fold 12 hpi) were elevated also they dominated at 6, 12, and 24 hpi (Fig 4C). Collectively these data suggest that, the significant induction of chemokines preferentially at early stage could be an indication of characteristic viral infection. The prominent regulation of the apoptosis related genes together with inflammatory chemokines at early stage of infection may indicate that many important events related to the virus pathogenicity had happened during early phase of infection. Understanding the cell signaling profiles driven by DE genes can be an indication of interaction between virus and host-cell as well as prediction of cell fate.

Selected DE genes form diverse gene ontologies showed consistent trend in RT-qPCR assay

To validate the data analyzed using microarray, expression of eight DE genes from diverse functional groups including chemokine (CXCL8, CCL20, CXCL3), anti-apoptotic (BCL2A1), transcription factor (FOSL1, JUN, and EGR1) and signal transduction (DUSP6 and GNB4) was tested by RT-qPCR using same RNA samples used for microarray hybridization at 6, 12, 24 and 48 hpi. When both RT-qPCR and microarray data showed a similar direction in fold change, the genes were considered validated. Fig 5 displays the detailed analysis of 8 DE genes by RT-qPCR. In consistent with microarray, EGR1 gene showed highest fold induction based on RT-qPCR analysis (69.8 fold at 12 hpi). The results indicated that the induction of CXCL8, CCl20, CXCL3, BCL2A1, FOSL1, JUN, EGR1 and DUSP6 showed consistent fold change direction both in RT-qPCR and microarray analysis except for GNB4. However, we did not find any significant differences in the level of expression of IFN-β and IFN-λ in both microarray as well as RT-qPCR analyses.

Discussion

The current study provides comprehensive gene expression analysis of HRV infection based on four time points using an in vitro system. We hypothesized that a discrete gene subsets could be disrupted along with development of infection. Undifferentiated epithelial cell culture model was selected in order to increase susceptibility of the cells to virus in order to reach uniform infection. The advantage of this system allows analysis of data mainly among infected cells rather than uninfected cells. Although replication of the rhinovirus can be low grade with non-cytopathic effect, in the current study HRV72 infection was determined by progression of characteristic CPE which was accompanied by efficient intracellular (cell lysate) and extracellular (cell culture supernatant) accumulation of the viral genome. This finding further provides an in vitro evidence for implication of HRV-B in LRTI. The transcriptional response to HRV infection was characterized by significant up- and down-regulation of the genes. The prominent down-regulation of the genes and low induction of antiviral genes could be related to the cleavage of transcriptional factors observed in Picornaviruses.

We aimed to identify the genes induced at early infection stage since modifications at early time points is mainly due to direct induction through viral signaling rather than feedback action of the other cell products which may play important role on severity and outcome of disease. In contrast to previous gene expression studies, our microarray analysis showed a significant difference between HRV-infected and mock-infected cells at 6 hpi. Gene expression analysis of air-liquid differentiated primary human tracheobronchial epithelial cells infected with HRV major and minor groups showed no DE genes at early infection stage. Similar result was observed by Proud et al. using in vivo system. This discrepancy might be due to patchy pattern of HRV infection in vivo and low susceptibility of differentiated epithelial cells for HRV infection.

The transcriptional response to HRV infection was further characterized by up-regulation of chemotactic factors at early stage. The monophasic nature of the chemokine induction started at 6 hpi and continued up to 12 hpi, potentially rejects the presence of secondary extracellular stimulus. The potential role of respiratory viruses for induction of airway inflammation has been demonstrated for HRV as well as other viruses. The induction of chemotactic factors in this study suggests a potential role of the HRV72 infected alveolar epithelial cells to initiate pulmonary inflammation, the possible role that has been demonstrated for primary human tracheobronchial epithelial cells infected with HRV-A strains. In addition, strain-specific differences has been reported among influenza viruses regarding the infectivity grade and inflammatory response. The current study showed that HRV72 induced expression of CXC (ELR-CXC) and CC chemokines and further suggests that neutrophilic and non-neutrophilic arms of the immune response could be activated upon HRV72 infection. Indeed uncontrolled expression of the chemokines can eventually lead to pathological inflammatory reactions.

The association between CXCL8 and reduced lung function, severity of asthma and airway hyper-responsiveness and allergic inflammation has been demonstrated in previous studies. Our findings further support the notion that HRV-induced release of CXCL8 and other members of CXC chemokine family (CXCL 1, 3 and 5) from alveolar epithelial cells could function as neutrophil chemoatracting signal and that perpetuate the pathogenesis of the HRV-B in LRTIs. CC chemokines including CCL2 and CCL20 were significantly induced upon HRV72 infection suggesting that non-neutrophilic immune response could be potentially provoked by this virus. Association of CCL2 level with respiratory symptoms has been shown in allergic individuals. Previous studies demonstrated that CCL2 was implicated in HRV induced airway hyper-responsiveness as a mechanism underlying asthma exacerbation. Further in vivo studies would address the relationship between immune parameters e.g. chemokine level and severity of the disease.

Contrary to expectations that induction of type I IFN and subsequent expression of interferon-stimulated genes (ISGs) are common phenomenon in HRV infected primary human tracheobronchial epithelial cells, this study did not find a significant induction of IFN mRNA expression. The induction of IFN-β mRNA and activation of Jak/STAT pathway was not also reported for HRV 14. It has been shown that IRF3 acts as an intermediate molecule to stimulate the induction of IFN-I and variety of ISGs in airway epithelial cells infected with HRV group A. Low level of the IFN-I in the study by Kotla et l. (2008) could be explained by the impaired activation of the IRF3 in the presence of the activated NF-κB and ATF-2. Variable grades of IFN type I induction in the strain-specific differences manner have been shown among influenza viruses and host antiviral which could play a critical role in clinical disease. The discrepancy in the induction or inhibition of IFN-I observed in HRV infections is not clear and could be different between primary and immortalized cells as well as different virus serotypes. This finding further emphasizes the important implication of developing IFN based treatment for HRV-B infections. Defective induction of major antiviral genes could be primarily related to the defective induction of IFN genes. HRV72 induced expression of tripartite motif containing (TRIM 10, 68, 26, and 69) genes, which play a broad role in innate immune response. The positive regulatory effect versus negative regulatory effect of TRIM proteins on antiviral response has been reported in other viruses and it is noteworthy to investigate the effect of TRIMs found in this study on HRV infection.

Viruses have evolved strategies to exploit cell signaling pathways for virus replication. Targeting of these pathways might be a valuable antiviral approach. MAPK signaling pathways play critical role in innate and adaptive immune response. Association of MAPK activation to the production of chemokines and cytokines has been shown in influenza infection[42]. Inactivation of MAPK pathways through DUSP also plays an important role in immune system[43]. Interestingly, DUSP6 was induced at first three time points in our study. Understanding the role of DUSP6 in MAPKs regulation in the context of HRV infection may require further investigations. Activation of signal transduction pathways could result in activation of transcription factors such as AP-1 and NF-κB. HRV72 induced expression of FOSL1, FOSL2, FOS and JUN genes. Combinations of the Fos and Jun proteins constitute dimeric transcription factor called AP-1. The AP-1 mediated production of chemokies has been demonstrated by Hoffmann et al. who delineated different roles for Fra-1 and c-Fos in transcription of IL8. Cooperative role of AP-1 in conjugation with NF-κB in IL8 induction has been shown for RSV[45]. More research on this topic including protein level expression as well as determination of absolute or cooperative interaction of AP-1 with other transcription factors in stimulation of chemokines could be useful for the develop of combinatorial therapy for HRV infection. Surprisingly EGR-1 was the DE gene expressed at the highest level in this study. Previous study reveal that HSV-1 viral load and subsequent mortality is reduced in EGR-1 knockout mice suggestive of EGR-1 role in regulation of virus replication. Further investigations would address the association between EGR-1 and HRV72 replication.

Genes involved in regulation of pro- versus anti-apoptotic responses were differentially expressed upon HRV infection in airway epithelial cells. The previous studies of apoptotic pathways indicated that different HRV strains can activate both intrinsic and extrinsic apoptotic pathways through caspase 9 and caspase 8, respectively. BCL2A and MCL1 belonged to the anti-apoptotic regulators of Bcl-2 family were significantly induced in response to HRV infection. The cytoprotective effect of over-expressed BCL2A and MCL1 genes played critical role in neutrophil survival and inhibited oxidant-induced necrosis and apoptosis in lung epithelial cells[51]. The induction of BCL2A1 and MCL-1 may be causing higher cell viability despite incomplete preservation of cell morphology and development of CPE in the current study. Over-expressed Bcl-2 and Bcl-xl in Coxsackievirus B3 infected cells preserved high cell viability and inhibited release of Cyt c from mitochondria. From a translational point of view, in depth understanding of pro/anti-apoptotic pathways regulated by the virus could lead to the development of novel therapeutics against HRV.

The category of the genes related to the wound healing and coagulation cascade were differentially expressed between HRV infected vs. un-infected airway epithelia cells. The expression of fibrinogen chains, α, β and γ, was down-regulated during the infection course, although the finding was different from previous studies using primary human epithelial cells. The over expression of fibrinogen was demonstrated in infected lung epithelial cells which could contribute to pathology of lung damage. PLAU/PLAUR system play complex role in airway system such as cellular adhesion, proliferation, migration and tissue remodeling. PLAU/PLAUR system was induced in current study, the finding that is inconsistent with previous reports. PLAUR plays a role in infiltration of immune cell is lung bacterial infection although didn’t show similar role in RSV and influenza using murine models. Coagulation plays an essential role in control and clearance of the infection. The inta-alveolar deposition of fibrin by epithelial cells is a typical response to the lung injury and infection[57]. On the other hand, the expression of the coagulation cascade related genes including THBD and MMP were induced here.

The study showed detail analysis of HRV-B/cell interaction using A549 cells derived from a human alveolar carcinoma cells with properties of type II alveolar epithelial cells. The finding based on the model, however, must be carefully interpreted and cautiously extrapolated to the real airway system in human. It is not clear to what extent HRV-B induced inflammatory responses are involved in protective or immunopathogenic features of the infection. The relative role of the immune factors in clearing versus pathology of the virus may also have like to patient underlying condition. The chemokine role in viral infection can be further characterized by analysing chemokine knockout mice in the context of infection. It is intriguing to speculate that defective elimination of HRV-B72 in the current in vitro system might be presented as a severe case in in vivo system. Further research regarding to the polymorphism of chemokine genes would be advantageous in determining genotype-phenotype associations of HRV infections and complications.

In summary, comprehensive gene expression analysis of human lung epithelial cells infected with HRV-B was presented in current study. Various functional gene categories were enriched at different time points post infection, which might be specific for HRVB species. This study provides a foundation for the further investigation of the genes and pathways that might be involved in the pathogenesis of the HRVB and the possible targets for treatment.