key: cord-0002725-w77jb4s8
authors: Reza Etemadi, Mohammad; Ling, King-Hwa; Zainal Abidin, Shahidee; Chee, Hui-Yee; Sekawi, Zamberi
title: Gene expression patterns induced at different stages of rhinovirus infection in human alveolar epithelial cells
date: 2017-05-30
journal: PLoS One
DOI: 10.1371/journal.pone.0176947
sha: 91b3ef8e5827c66d841a773d4d1190b6007f4735
doc_id: 2725
cord_uid: w77jb4s8

Human rhinovirus (HRV) is the common virus that causes acute respiratory infection (ARI) and is frequently associated with lower respiratory tract infections (LRTIs). We aimed to investigate whether HRV infection induces a specific gene expression pattern in airway epithelial cells. Alveolar epithelial cell monolayers were infected with HRV species B (HRV-B). RNA was extracted from both supernatants and infected monolayer cells at 6, 12, 24 and 48 hours post infection (hpi) and transcriptional profile was analyzed using Affymetrix GeneChip and the results were subsequently validated using quantitative Real-time PCR method. HRV-B infects alveolar epithelial cells which supports implication of the virus with LRTIs. In total 991 genes were found differentially expressed during the course of infection. Of these, 459 genes were up-regulated whereas 532 genes were down-regulated. Differential gene expression at 6 hpi (187 genes up-regulated vs. 156 down-regulated) were significantly represented by gene ontologies related to the chemokines and inflammatory molecules indicating characteristic of viral infection. The 75 up-regulated genes surpassed the down-regulated genes (35) at 12 hpi and their enriched ontologies fell into discrete functional entities such as regulation of apoptosis, anti-apoptosis, and wound healing. At later time points of 24 and 48 hpi, predominated down-regulated genes were enriched for extracellular matrix proteins and airway remodeling events. Our data provides a comprehensive image of host response to HRV infection. The study suggests the underlying molecular regulatory networks genes which might be involved in pathogenicity of the HRV-B and potential targets for further validations and development of effective treatment.

Human rhinovirus (HRV), a non-segmented positive sense RNA from Picornaviridae, [1] , is the most prevalent human respiratory pathogen [2] , consisted of three species, HRV-A, B and C [3] . HRV is associated with common colds and concomitant sinusitis [4] and otitis media episodes [5] as well as lower respiratory tract infection (LRTI) [6] . 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 [7, 8] . Association of HRV with LRTIs has also been demonstrated in lower airway samples [9, 10] . Human airway epithelial cell is the primary target for HRV infection leading to interstitial pneumonia [11, 12] . 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 [13, 14] . Global transcriptional response of HRV-infected pulmonary epithelial cells has been investigated using in vitro [15, 16] and in vivo [17] 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 [18, 19] , the transcriptomic data is rather limited and have not been extensively understood with respect to different species [20] . 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.

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 CO 2 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.

H1Hela cells were plated over 6-well plates at a cell density of 4-5 × 10 5 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 [21, 22] . The titer of the purified virus was 5.2 × 10 9 PFU/mL, as depicted in S2 Fig.

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 [23] . A549 cells where seeded in 25 cm 2 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.

The A549 cells grown in 25 cm 2 culture flasks were infected with HRV72 at multiplicity of infection (MOI) of 10 or mock-infected at three independent replicates with a week interval [24] . 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.

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.

Three independent replicates of RNA samples from mock-infected and HRV-infected A549cells 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 GeneChipPrime-View 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).

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/microarrayanalysis/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 [25] 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 0 /5 0 ratios of internal control probe sets including GAPDH and ACTB were less than 3 for all the arrays indicating appropriate sample quality [26] . 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 log 2 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 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.

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.

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. Oneway 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.

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 [15] . 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. 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 [19] .

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 

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 ( 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.

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.

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 [27] . 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 noncytopathic effect [28] , 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 [11, 28] . 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 [29] . 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 [15, 17] , 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 [15] . Similar result was observed by Proud et al. using in vivo system [17] . 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 [24] . The potential role of respiratory viruses for induction of airway inflammation has been demonstrated for HRV [16] as well as other viruses [30] . 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 [15, 17] . In addition, strain-specific differences has been reported among influenza viruses regarding the infectivity grade and inflammatory response [31] . 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 [32] and allergic inflammation [33] has been demonstrated in previous studies. Our findings further support the notion that HRV-induced release of CXCL8 [32, 34] and other members of CXC chemokine family (CXCL 1, 3 and 5) [35] 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 [36] . Previous studies demonstrated that CCL2 was implicated in HRV induced airway hyper-responsiveness as a mechanism underlying asthma exacerbation [36, 37] . 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 [15, 17] , 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 [20, 38] . 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 [39] . Low level of the IFN-I in the study by Kotla et l. (2008) [38] 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 [31, 40] . 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 [41] .

Viruses have evolved strategies to exploit cell signaling pathways for virus replication. Targeting of these pathways might be a valuable antiviral approach [42] . 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. [44] 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 [46] . 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 [47] and caspase 8, respectively [48] . BCL2A and MCL1 belonged to the anti-apoptotic regulators of Bcl-2 family [49] were significantly induced in response to HRV infection. The cytoprotective effect of over-expressed BCL2A and MCL1 genes played critical role in neutrophil survival [50] 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 [52] . 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 [15, 17] . The over expression of fibrinogen was demonstrated in infected lung epithelial cells which could contribute to pathology of lung damage [53, 54] . PLAU/PLAUR system play complex role in airway system such as cellular adhesion, proliferation, migration and tissue remodeling [55] . PLAU/PLAUR system was induced in current study, the finding that is inconsistent with previous reports [15, 17] . 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 [56] . 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 [58] .

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.

Supporting information S1 Fig. SDS page analysis of purified HRV72 virus. 20 μl of each sucrose fractions were mixed with 20 μl of 1× SDS page sample buffer and denatured for 10 min at 100˚C followed by gel running and subsequent staining of SDS page with coomassie blue. Sucrose fractions from up to down were run onto the gel from right to left order. In the UP6 fraction, the purified virion shows three distinct bands which is more likely corresponded to vp1-3. (TIF) S2 Fig. HRV72 plaque formation on H1Hela cells. H1Hela cells were seeded in 6-well plates at a cell density of 4-5 × 10 5 cells per well and incubated overnight at 37˚C. Wells infected with 0.1 ml of dilutions of virus ranging from 10 1 (left up) to 10 −10 and a mock-infected (right down) in duplicates. At the end of 5-day incubation period at 34˚C, the carboxymethylcellulose over layer was decanted and stained with 0.1% crystal violet solution. Plaques were counted in the 10 −7 and 10 −8 wells using this equation: PFU/ml = (Mean of plaques / amount of inoculum (ml)) / well dilution. (TIF) S3 Fig. Toxicity of HRV72 on A549 cells. Cells were seeded at optimal density of 1×10 4 cells/ well in triplicates in 96 well plate and the cell infection was carried out at 24 h when cells were almost 80% confluent. The cells were inoculated with various MOI (1, 5, 7, 10, 15, 20, 40 and 100 pfu/mL) and incubated at 34˚C with 5% CO 2 . Toxicity was measured by determining metabolic activity of the A549 cells using MTT assay. Data are presented as average of three independent experiments. The cell viability was calculated with this equation: ¼ absorbance of sample absorbance of negative control  100. The absorbance was measured at 570 nm. (TIF) S1 

Picornaviridae: The Viruses and Their Replication

Role of respiratory viruses in acute upper and lower respiratory tract illness in the first year of life: a birth cohort study

Proposals for the classification of human rhinovirus species A, B and C into genotypically assigned types

Detection of rhinovirus in sinus brushings of patients with acute community-acquired sinusitis by reverse transcription-PCR

Virological and serological analysis of rhinovirus infections during the first two years of life in a cohort of children

Rhinovirus is an important pathogen in upper and lower respiratory tract infections in Mexican children

Do rhinoviruses cause pneumonia in children?

Rhinovirus infections in children: A retrospective and prospective hospital-based study

Exacerbations of Asthma in Adults during Experimental Rhinovirus Infection 1-3. American Review of Respiratory Disease

Detection of rhinovirus RNA in lower airway cells during experimentally induced infection

Pneumonia Caused by Rhinovirus

Similar Frequency of Rhinovirus-Infectible Cells in Upper and Lower Airway Epithelium

Neutrophil degranulation and cell lysis is associated with clinical severity in virus-induced asthma

Role of rhinovirus infections in asthma

Rhinovirus induces airway epithelial gene expression through double-stranded RNA and IFN-dependent pathways

Rhinovirusinduced modulation of gene expression in bronchial epithelial cells from subjects with asthma

Gene Expression Profiles during In Vivo Human Rhinovirus Infection

Clinical and molecular epidemiology of human rhinovirus C in children and adults in Hong Kong reveals a possible distinct human rhinovirus C subgroup

Biodiversity and clinico-demographic characteristics of human rhinoviruses from hospitalized children with acute lower respiratory tract infections in Malaysia

Human Rhinovirus Attenuates the Type I Interferon Response by Disrupting Activation of Interferon Regulatory Factor 3

Major and minor group rhinoviruses elicit differential signaling and cytokine responses as a function of receptor-mediated signal transduction

Rhinovirus plaque formation in WI-38 cells with methylcellulose overlay

Understanding cell toxicology: principles and practice Ellis Horwood

Role of p38 Mitogen-Activated Protein Kinase in Rhinovirus-Induced Cytokine Production by Bronchial Epithelial Cells

Summaries of Affymetrix Gene-Chip probe level data

Genechip Expression Analysis Technical Manual

Resistance of differentiated human airway epithelium to infection by rhinovirus

Low Grade Rhinovirus Infection Induces a Prolonged Release of IL-8 in Pulmonary Epithelium

Cytopathogenesis and inhibition of host gene expression by RNA viruses

Expression of respiratory syncytial virus-induced chemokine gene networks in lower airway epithelial cells revealed by cDNA microarrays

Infection of lung epithelial cells with pandemic 2009 A (H1N1) influenza viruses reveals isolate-specific differences in infectivity and host cellular responses

Effect of experimental rhinovirus 16 colds on airway hyperresponsiveness to histamine and interleukin-8 in nasal lavage in asthmatic subjects in vivo

Interleukin-8 is a potent eosinophil chemo-attractant

Rhinovirus-infected epithelial cells produce more IL-8 and RANTES compared with other respiratory viruses

Constitutive and stimulated MCP-1, GRO alpha, beta, and gamma expression in human airway epithelium and bronchoalveolar macrophages

Nasal cytokine responses to natural colds in asthmatic children

Macrophage/ epithelial cell CCL2 contributes to rhinovirus-induced hyperresponsiveness and inflammation in a mouse model of allergic airways disease

Attenuation of the type I interferon response in cells infected with human rhinovirus

Role of doublestranded RNA pattern recognition receptors in rhinovirus-induced airway epithelial cell responses

Variation in the ability of human influenza A viruses to induce and inhibit the IFN-β pathway

Toll-like receptors and their crosstalk with other innate receptors in infection and immunity

Influenza viruses and MAP kinase cascades-Novel targets for an antiviral intervention?

DUSP meet immunology: dual specificity MAPK phosphatases in control of the inflammatory response

MEK1-dependent Delayed Expression of Fos-related Antigen-1 Counteracts c-Fos and p65 NF-κB-mediated Interleukin-8 Transcription in Response to Cytokines or Growth Factors

Activator protein-1 is the preferred transcription factor for cooperative interaction with nuclear factor-κB in respiratory syncytial virus-induced interleukin-8 gene expression in airway epithelium

Suppression of transcription factor early growth response 1 reduces herpes simplex virus lethality in mice

Apoptotic events induced by human rhinovirus infection

NuMA and nuclear lamins are cleaved during viral infection-inhibition of caspase activity prevents cleavage and rescues HeLa cells from measles virus-induced but not from rhinovirus 1B-induced cell death

Apoptosis: a review of programmed cell death

BCL-2 family expression in human neutrophils during delayed and accelerated apoptosis

Bcl-2-related protein A1 is an endogenous and cytokine-stimulated mediator of cytoprotection in hyperoxic acute lung injury

Bcl-2 and Bcl-xL overexpression inhibits cytochrome c release, activation of multiple caspases, and virus release following coxsackievirus B3 infection

The severe acute respiratory syndrome coronavirus 3a protein up-regulates expression of fibrinogen in lung epithelial cells

Extrahepatic expression of the γ chain of fibrinogen during an acute phase response

uPAR: a versatile signalling orchestrator

Urokinase receptor-deficient mice mount an innate immune response to and clarify respiratory viruses as efficiently as wild-type mice

Novel role of the human alveolar epithelium in regulating intra-alveolar coagulation

The interleukin 6− 174 C/C genotype predicts greater rhinovirus illness

The authors thank Research Instruments Sdn. Bhd., Malaysia, for providing facilities and technical support in performing microarray test. We would like to extent our appreciation to the members of Medical Genetics Laboratory, Department of Biomedical Sciences for allowing us to carry out part of above study. We thank Sara Ansari, from Department of Anatomy, for excellent technical assistance.