# A Cross-Study Biomarker Signature of Human Bronchial Epithelial Cells Infected with Respiratory Syncytial Virus

Authors: Luiz Gustavo Gardinassi
Journal: Advances in Virology
 (2016)

---

## Abstract

Respiratory syncytial virus (RSV) is a major cause of lower respiratory tract infections in children, elderly, and immunocompromised individuals. Despite of advances in diagnosis and treatment, biomarkers of RSV infection are still unclear. To understand the host response and propose signatures of RSV infection, previous studies evaluated the transcriptional profile of the human bronchial epithelial cell line—BEAS-2B—infected with different strains of this virus. However, the evolution of statistical methods and functional analysis together with the large amount of expression data provide opportunities to uncover novel biomarkers of inflammation and infections. In view of those facts publicly available microarray datasets from RSV-infected BEAS-2B cells were analyzed with linear model-based statistics and the platform for functional analysis InnateDB. The results from those analyses argue for the reevaluation of previously reported transcription patterns and biological pathways in BEAS-2B cell lines infected with RSV. Importantly, this study revealed a biosignature constituted by genes such asABCC4,ARMC8,BCLAF1,EZH1,FAM118A,FAM208B,FUS,HSPH1,KAZN,MAP3K2,N6AMT1,PRMT2,S100PBP,SERPINA1,TLK2,ZNF322, andZNF337 which should be considered in the development of new molecular diagnosis tools.


---

## Body

## 1. Introduction

Respiratory syncytial virus (RSV) is a major etiologic agent causing acute lower respiratory infections that can progress to bronchiolitis and pneumonia in children, elderly, and immunocompromised individuals [1, 2]. RSV outbreaks are influenced by virus diversity and evolution [3, 4], environmental factors [5], and host immunity [6].The epithelium is the primary site for host-virus interface, where cells recognize pathogen-associated patterns on microbes through innate immunity receptors [7, 8]. Indeed, epithelial cells constitute an important line of defense against RSV and other airborne pathogens [9]. They form a physical barrier and produce mucus to inhibit microbes from entering the body. Moreover, they express molecules with antimicrobial properties, as lysozyme, lactoferrin, collectins, and antimicrobial peptides [10]. Two human cell lines have been extensively used to understand the interaction between host and RSV, the alveolar epithelial cell, A549, and one from proximal airways, the bronchial epithelial cell, BEAS-2B.Genome-wide microarrays are powerful tools to investigate host transcriptional response during infections in the pulmonary epithelium, including those induced by RSV [11, 12]. Indeed, two studies evaluated the patterns of gene expression from BEAS-2B cell lines infected with RSV [10, 13]. However, it is intriguing that after 4 h of infection Huang and collaborators (2008) found that RSV-modulated genes were only associated with the neuroactive ligand-receptor interaction pathway [13]; in contrast, Mayer and collaborators (2007) identified that the same time of RSV infection of BEAS-2B cells induced transcriptional changes similar to those found for other respiratory pathogens asPseudomonas aeruginosa [10]. In spite of differences, publicly available microarray data offers an interesting opportunity to reveal common features of RSV induced transcriptional profiles to understand the early response of BEAS-2B cell lines and extend the knowledge on biomarkers of acute infections with this virus. Therefore, those datasets were evaluated in a meta-analysis by fitting linear models for each array probe and Empirical Bayesian approach to detect transcriptional changes that revealed significant associations with unreported pathways. Of importance, this strategy also rendered a biomarker signature of BEAS-2B cell lines infected with RSV that can be useful for the design of molecular diagnosis tools.


## 2. Materials and Methods

The datasets GSE3397 and GSE6802 were obtained from GEO database (http://www.ncbi.nlm.nih.gov/), which compared BEAS-2B cells infected with RSV with control experiments. Only arrays in which cells were infected with RSV for 4 h were selected for further analysis. Raw data were processed using the R Language and Environment for Statistical Computing (R) 3.2.0 [14] and Bioconductor 3.1 [15]. Theaffy package for R [16] was used to perform quality control when applicable. Data was log
2
⁡ transformed and quantile normalization was applied for dataset GSE3397 due the absence of CEL files. The dataset GSE6802 was already RMA normalized. Batch effects were corrected with Combat( ) function [17] ofsva package for R [18]. Expression data were weighted with the arrayWeights( ) function fromlimma package for R [19]. Differential gene expression was also evaluated withlimma package for R [19], whereby differentially expressed genes (DEGs) were identified by a false discovery rate (FDR) <0.05. Hierarchical clustering was performed with Euclidian distance for metric calculations and the complete linkage method, which were displayed as heatmaps drawn withgplots package for R [20]. Pathway analyses were performed with the online platform for functional analysis InnateDB [21] and significant pathway overrepresentation was computed with hypergeometrical distribution and Benjamini-Hochberg correction for multiple comparisons. Significantly enriched pathways were determined by a P value < 0.05 and FDR < 0.1.


## 3. Results and Discussion

### 3.1. Dataset Selection and Preprocessing Analysis

To define a robust transcriptional signature of BEAS-2B acutely infected with RSV, two publicly available datasets, GSE3397 and GSE6802, were used to conduct a meta-analysis from which data were extracted for BEAS-2B cells infected with RSV for 4 h and controls. First, background subtracted expression data from GSE3397 (Figure1(a)) were preprocessed and normalized (Figure 1(b)). However, in a first attempt to conduct differential gene expression analysis usinglimma [19], there were no statistically significant differences in gene expression. Therefore, principal component analysis (PCA) was used to evaluate the expression profiles of each array and, except for arrays named here Control2 and RSV2, the consistent pattern of clustering in Figure 1(c) suggests a batch effect. After normalization, this effect was even more evident (Figure 1(d)), which led to the speculation that Huang and collaborators (2008) [13] analyzed only three microarray experiments from this dataset based on the assumption that differences found for those microarrays were due to failures in experimental procedures; however they did not consider or correct for batch effects. In view of those facts, the datasets were adjusted with Combat function for R, which removed such effects from GSE3397 expression data (Figure 1(e)). Batch correction of GSE3397 did not change the profiles of arrays Control2 and RSV2; nevertheless, those arrays were included in further analysis because the variation observed in this experiment could have a substantial impact over the final result. Even adverse experimental variations that may change the overall expression patterns of a dataset could be useful to power up the identification of genes that are robustly modulated in BEAS-2B cells infected with RSV. The expression dataset GSE6802 (Figure 1(f)) was also included in the analysis. PCA from expression data extracted from GEO demonstrates that most of the variability between the arrays is explained (76.6%) by the infection with RSV, as the standardized PC1 separates RSV-infected from control arrays (Figure 1(g)), whereas standardized PC2 (11.4%) separates one pair of arrays (RSV_3 and ctrl2) and, although these arrays are supposedly from different batches, clustering features of this axis also suggested a batch effect (Figure 1(g)). log2 transformation of data impacted the profile of array RSV_1 however did not change the profiles from RSV_3 and ctrl_2 (Figure 1(h)). Combat( ) function was also applied to the expression dataset GSE6802; however, PCA shows that the adjustment did not to improve further clustering between specific arrays (Supplementary Figure  1; see Supplementary Material available online at http://dx.doi.org/10.1155/2016/3605302). In view of that, downstream analyses were carried out with normalized log2 transformed data.Figure 1
Preprocessing analysis of GEO datasets GSE3397 and GSE6802. (a) Boxplot of GSE3397, raw expression data. (b) Boxplot of GSE3397, normalized expression data. (c) Principal component analysis of GSE3397, raw expression data. (d) Principal component analysis of GSE3397, normalized expression data. (e) Principal component analysis of GSE3397, normalized and batch corrected expression data. (f) Boxplot of GSE6802, RMA normalized expression data. (g) Principal component analysis of GSE6802, normalized expression data. (h) Principal component analysis of GSE6802, log2 transformed RMA normalized expression data.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


### 3.2. Differential Gene Expression

Next, linear model-based statistical analyses with a FDR < 0.05 were conducted to identify differentially expressed genes (DEGs). The dataset GSE3397 exhibited ninety-four DEGs (Figure2(a) and Table 1). Those genes are highly discordant from DEGs previously reported by Huang and collaborators (2008) [13], which identified 277 DEGs based on different statistical analysis and assumptions. Fifty genes were downregulated and forty-four were upregulated (Table 1). The differences found in this study might reflect the inclusion of all microarray experiments from controls and 4 h after RSV infection; exclusion of expression data from 24 h after RSV infection; distinct preprocessing approaches as normalizing method and batch effect correction; and the assessment of statistical significance with a linear model-based method and corrected P values. In contrast, 1965 DEGs were identified for the dataset GSE6802. The top hundred DEGs ranked by fold changes (Figure 2(b) and Table 2) included genes such asJUNB,KLF4,CXCL1,CXCL2, andIL6, which are in agreement with those reported by Mayer and collaborators (2007) [10]. Several factors should account for the notable differences in expression analysis from both datasets. First, different RSV strains were used to stimulate BEAS-2B cells. Second, experimental conditions of controls were also different, as control experiments from GSE3397 were incubated with vehicle (not specified) and those from GSE6802 were not stimulated. Third, despite both datasets being generated with affymetrix microarray platform, those include distinct versions, HU133 plus 2.0 for GSE3397 and HU133A 2.0 for GSE6802.Table 1
Differentially expressed genes identified in dataset GSE3397.
ProbeID
Gene symbol
Gene name
l
o
g
2 fold change
FDR
1560754_at
CMTM7
CKLF like MARVEL transmembrane domain containing 7
−1,54756
0,017104
239439_at
AFF4
AF4/FMR2 family member 4
−1,53581
0,023832
238929_at
SRSF8
Serine/arginine-rich splicing factor 8
−1,51887
0,018433
223142_s_at
UCK1
Uridine-cytidine kinase 1
−1,47939
0,017104
242636_at
PRCP
Prolylcarboxypeptidase
−1,45095
0,034358
228007_at
CEP85L
Centrosomal protein 85 kDa-like
−1,4103
0,017104
235573_at
HSPH1
Heat shock protein family H (Hsp110) member 1
−1,39959
0,0371
228391_at
CYP4V2
Cytochrome P450 family 4 subfamily V member 2
−1,38799
0,01671
219376_at
ZNF322
Zinc finger protein 322
−1,3491
0,046761
1553689_s_at
METTL6
Methyltransferase like 6
−1,34723
0,017104
242837_at
SRSF4
Serine/arginine-rich splicing factor 4
−1,34071
0,044693
237215_s_at
TFRC
Transferrin receptor
−1,32685
0,017104
208819_at
RAB8A
RAB8A, member RAS oncogene family
−1,32593
0,042264
236665_at
CCDC18
Coiled-coil domain containing 18
−1,31494
0,034201
206147_x_at
SCML2
Sex comb on midleg-like 2 (Drosophila)
−1,30586
0,016454
229325_at
ZZZ3
Zinc finger ZZ-type containing 3
−1,30495
0,017104
1565716_at
FUS
FUS RNA binding protein
−1,29415
0,049505
205062_x_at
ARID4A
AT-rich interaction domain 4A
−1,28877
0,033039
1552312_a_at
MFAP3
Microfibrillar associated protein 3
−1,28521
0,046511
223223_at
ARV1
ARV1 homolog, fatty acid homeostasis modulator
−1,27987
0,023832
232001_at
PRKCQ-AS1
PRKCQ antisense RNA 1
−1,27987
0,035983
233195_at
DNAI1
Dynein axonemal intermediate chain 1
−1,25963
0,047083
219094_at
ARMC8
Armadillo repeat containing 8
−1,25527
0,043392
235232_at
GMEB1
Glucocorticoid modulatory element binding protein 1
−1,2492
0,046511
218643_s_at
CRIPT
CXXC repeat containing interactor of PDZ3 domain
−1,24229
0,0371
1566851_at
TRIM42
Tripartite motif containing 42
−1,24057
0,042149
221821_s_at
KANSL2
KAT8 regulatory NSL complex subunit 2
−1,23799
0,017104
244115_at
FAM126A
Family with sequence similarity 126 member A
−1,23114
0,033039
215541_s_at
DIAPH1
Diaphanous related formin 1
−1,22774
0,033039
203196_at
ABCC4
ATP binding cassette subfamily C member 4
−1,22519
0,033039
225024_at
RPRD1B
Regulation of nuclear pre-mRNA domain containing 1B
−1,22264
0,043765
37860_at
ZNF337
Zinc finger protein 337
−1,22095
0,023832
212997_s_at
TLK2
Tousled like kinase 2
−1,21841
0,04814
225690_at
CDK12
Cyclin-dependent kinase 12
−1,21083
0,0371
232103_at
BPNT1
3′(2′), 5′-Bisphosphate nucleotidase 1
−1,20748
0,0371
224848_at
CDK6
Cyclin-dependent kinase 6
−1,20247
0,0371
214962_s_at
NUP160
Nucleoporin 160 kDa
−1,20247
0,046319
219629_at
FAM118A
Family with sequence similarity 118 member A
−1,19831
0,028374
212290_at
SLC7A1
Solute carrier family 7 member 1
−1,19748
0,042264
227187_at
CBLL1
Cbl proto-oncogene like 1, E3 ubiquitin protein ligase
−1,19582
0,030047
233208_x_at
CPSF2
Cleavage and polyadenylation specific factor 2
−1,19334
0,046319
230566_at
MORC2-AS1
MORC2 antisense RNA 1
−1,17691
0,0371
238795_at
FAM208B
Family with sequence similarity 208 member B
−1,17609
0,0371
204980_at
CLOCK
Clock circadian regulator
−1,17283
0,0371
238653_at
LRIG2
Leucine-rich repeats and immunoglobulin like domains 2
−1,17202
0,048527
229939_at
ENDOV
Endonuclease V
−1,16878
0,041349
218185_s_at
ARMC1
Armadillo repeat containing 1
−1,16151
0,046319
201083_s_at
BCLAF1
BCL2 associated transcription factor 1
−1,15509
0,049505
227840_at
C2orf76
Chromosome 2 open reading frame 76
−1,15109
0,042264
201686_x_at
API5
Apoptosis inhibitor 5
−1,14076
0,046761
221699_s_at
DDX50
DEAD-box helicase 50
1,140764
0,046511
1556178_x_at
TAF8
TATA-box binding protein associated factor 8
1,159096
0,034358
205623_at
ALDH3A1
Aldehyde dehydrogenase 3 family member A1
1,163927
0,049505
212495_at
KDM4B
Lysine demethylase 4B
1,193336
0,044693
1569057_s_at
MIA3
Melanoma inhibitory activity family member 3
1,193336
0,047866
222494_at
FOXN3
Forkhead box N3
1,19582
0,048527
223311_s_at
MTA3
Metastasis associated 1 family member 3
1,19582
0,041439
215424_s_at
SNW1
SNW domain containing 1
1,196649
0,049505
213478_at
KAZN
Kazrin, periplakin interacting protein
1,19914
0,025143
227864_s_at
MVB12A
Multivesicular body subunit 12A
1,201636
0,030287
228674_s_at
EML4
Echinoderm microtubule associated protein like 4
1,204137
0,040345
224196_x_at
DPH5
Diphthamide biosynthesis 5
1,205808
0,025143
224652_at
CCNY
Cyclin Y
1,207481
0,046761
212968_at
RFNG
RFNG O-fucosylpeptide 3-beta-N-acetylglucosaminyltransferase
1,211673
0,0371
1555486_a_at
PRR5L
Proline rich 5 like
1,212513
0,017104
232837_at
KIF13A
Kinesin family member 13A
1,214195
0,042264
224320_s_at
MCM8
Minichromosome maintenance 8 homologous recombination repair factor
1,217566
0,033039
230131_x_at
ARSD
Arylsulfatase D
1,221793
0,0371
218225_at
ECSIT
ECSIT signalling integrator
1,224336
0,034358
222610_s_at
S100PBP
S100P binding protein
1,226885
0,030047
32259_at
EZH1
Enhancer of zeste 1 polycomb repressive complex 2 subunit
1,229439
0,0371
203854_at
CFI
Complement factor I
1,232852
0,042264
221600_s_at
AAMDC
Adipogenesis associated, Mth938 domain containing
1,260503
0,0371
209558_s_at
HIP1R
Huntingtin interacting protein 1 related
1,263127
0,042264
224814_at
DPP7
Dipeptidyl peptidase 7
1,26488
0,016454
232280_at
SLC25A29
Solute carrier family 25 member 29
1,277214
0,030047
228424_at
NAALADL1
N-Acetylated alpha-linked acidic dipeptidase-like 1
1,286989
0,042264
203409_at
DDB2
Damage specific DNA binding protein 2
1,288775
0,023832
229975_at
BMPR1B
Bone morphogenetic protein receptor type 1B
1,297739
0,034358
227073_at
MAP3K2
Mitogen-activated protein kinase kinase kinase 2
1,297739
0,017104
225347_at
ARL8A
ADP ribosylation factor like GTPase 8A
1,298639
0,02672
221774_x_at
SUPT20H
SPT20 homolog, SAGA complex component
1,308578
0,016454
223679_at
CTNNB1
Catenin beta 1
1,318594
0,018043
227679_at
HDAC11
Histone deacetylase 11
1,328686
0,044693
220020_at
XPNPEP3
X-Prolyl aminopeptidase 3, mitochondrial
1,342573
0,031097
203199_s_at
MTRR
5-Methyltetrahydrofolate-homocysteine methyltransferase reductase
1,360371
0,017104
228722_at
PRMT2
Protein arginine methyltransferase 2
1,370783
0,016454
228951_at
SLC38A7
Solute carrier family 38 member 7
1,431969
0,016454
217529_at
ORAI2
ORAI calcium release-activated calcium modulator 2
1,453973
0,043775
220311_at
N6AMT1
N-6 adenine-specific DNA methyltransferase 1 (putative)
1,460032
0,017104
213402_at
ZNF787
Zinc finger protein 787
1,469169
0,017104
226055_at
ARRDC2
Arrestin domain containing 2
1,477338
0,017104
219756_s_at
POF1B
Premature ovarian failure, 1B
1,580083
0,016454
202833_s_at
SERPINA1
Serpin peptidase inhibitor, clade A (alpha-1 antiproteinase, antitrypsin), and member 1
2,488023
0,0371Table 2
Top hundred differentially expressed genes identified in dataset GSE6802.
ProbeID
Gene symbol
Gene name
l
o
g
2 fold change
FDR
212615_at
CHD9
Chromodomain helicase DNA binding protein 9
−3,69609
0,00131
221840_at
PTPRE
Protein tyrosine phosphatase, receptor type E
−3,56524
0,000195
220817_at
TRPC4
Transient receptor potential cation channel subfamily C member 4
−3,39168
0,001582
221703_at
BRIP1
BRCA1 interacting protein C-terminal helicase 1
−2,88786
0,021463
207012_at
MMP16
Matrix metallopeptidase 16
−2,82647
0,000119
219494_at
RAD54B
RAD54 homolog B (S. cerevisiae)
−2,81279
0,000177
207034_s_at
GLI2
GLI family zinc finger 2
−2,79723
0,005157
203518_at
LYST
Lysosomal trafficking regulator
−2,75872
5,90E − 05
205282_at
LRP8
LDL receptor related protein 8
−2,7549
0,000311
214440_at
NAT1
N-Acetyltransferase 1 (arylamine N-acetyltransferase)
−2,68515
0,001777
219627_at
ZNF767P
Zinc finger family member 767, pseudogene
−2,67957
0,00024
218984_at
PUS7
Pseudouridylate synthase 7 (putative)
−2,67586
0,001308
206554_x_at
SETMAR
SET domain and mariner transposase fusion gene
−2,63536
0,002432
219779_at
ZFHX4
Zinc finger homeobox 4
−2,62624
0,001411
213103_at
STARD13
StAR related lipid transfer domain containing 13
−2,57219
0,002525
210138_at
RGS20
Regulator of G-protein signaling 20
−2,55974
0,000415
204291_at
ZNF518A
Zinc finger protein 518A
−2,54383
9,70E − 05
204651_at
NRF1
Nuclear respiratory factor 1
−2,49147
0,003659
205408_at
MLLT10
Myeloid/lymphoid or mixed-lineage leukemia; translocated to, 10
−2,48975
5,10E − 05
219581_at
TSEN2
tRNA splicing endonuclease subunit 2
−2,45377
0,001774
218242_s_at
SUV420H1
Lysine methyltransferase 5B
−2,44698
0,000754
203242_s_at
PDLIM5
PDZ and LIM domain 5
−2,43851
0,001699
203868_s_at
VCAM1
Vascular cell adhesion molecule 1
−2,43513
0,000761
220206_at
ZMYM1
Zinc finger MYM-type containing 1
−2,36362
0,008439
207616_s_at
TANK
TRAF family member associated NFKB activator
−2,34567
0,000424
218303_x_at
KRCC1
Lysine-rich coiled-coil 1
−2,34567
0,003187
218490_s_at
ZNF302
Zinc finger protein 302
−2,32785
0,001816
206876_at
SIM1
Single-minded family bHLH transcription factor 1
−2,32624
0,001681
219128_at
C2orf42
Chromosome 2 open reading frame 42
−2,28628
0,002926
212861_at
MFSD5
Major facilitator superfamily domain containing 5
−2,27048
0,000823
218653_at
SLC25A15
Solute carrier family 25 member 15
−2,25636
0,000562
206943_at
TGFBR1
Transforming growth factor beta receptor I
−2,24856
0,025349
201995_at
EXT1
Exostosin glycosyltransferase 1
−2,247
0,000421
221430_s_at
RNF146
Ring finger protein 146
−2,23457
0,001084
212286_at
ANKRD12
Ankyrin repeat domain 12
−2,2253
0,00029
219544_at
BORA
Bora, aurora kinase A activator
−2,21914
0,000333
210455_at
R3HCC1L
R3H domain and coiled-coil containing 1 like
−2,2176
0,0039
219459_at
POLR3B
Polymerase (RNA) III subunit B
−2,2176
0,000832
219078_at
GPATCH2
G-patch domain containing 2
−2,19923
0,000723
204547_at
RAB40B
RAB40B, member RAS oncogene family
−2,17648
0,001741
209760_at
KIAA0922
KIAA0922
−2,17347
0,001048
218791_s_at
KATNBL1
Katanin regulatory subunit B1 like 1
−2,17347
0,001187
205173_x_at
CD58
CD58 molecule
−2,17196
0,00022
204352_at
TRAF5
TNF receptor associated factor 5
−2,16895
0,002659
212441_at
KIAA0232
KIAA0232
−2,16595
0,006084
204236_at
FLI1
Fli-1 proto-oncogene, ETS transcription factor
−2,15397
0,005141
203072_at
MYO1E
Myosin IE
−2,15248
0,000154
219904_at
ZSCAN5A
Zinc finger and SCAN domain containing 5A
−2,14801
0,00144
219133_at
OXSM
3-Oxoacyl-ACP synthase, mitochondrial
−2,12285
0,002424
205798_at
IL7R
Interleukin 7 receptor
−2,11257
0,00506
205476_at
CCL20
C-C motif chemokine ligand 20
4,613942
9,50E − 05
213497_at
ABTB2
Ankyrin repeat and BTB domain containing 2
4,623547
1,40E − 05
219179_at
DACT1
Dishevelled-binding antagonist of beta-catenin 1
4,642816
9,00E − 06
219228_at
ZNF331
Zinc finger protein 331
4,723971
6,00E − 06
213139_at
SNAI2
Snail family zinc finger 2
4,76673
1,40E − 05
218177_at
CHMP1B
Charged multivesicular body protein 1B
4,806544
1,00E − 05
203304_at
BAMBI
BMP and activin membrane-bound inhibitor
4,826576
3,00E − 06
201631_s_at
IER3
Immediate early response 3
4,833271
3,00E − 06
218559_s_at
MAFB
v-maf avian musculoaponeurotic fibrosarcoma oncogene homolog B
4,870264
0,000468
220266_s_at
KLF4
Kruppel-like factor 4 (gut)
4,890561
0,00022
209211_at
KLF5
Kruppel-like factor 5 (intestinal)
4,924578
0,002036
209681_at
SLC19A2
Solute carrier family 19 member 2
4,927992
5,90E − 05
205266_at
LIF
Leukemia inhibitory factor
4,955395
2,20E − 05
204790_at
SMAD7
SMAD family member 7
5,073566
0,000283
221667_s_at
HSPB8
Heat shock protein family B (small) member 8
5,422657
2,90E − 05
212665_at
TIPARP
TCDD-inducible poly(ADP-ribose) polymerase
5,525098
1,00E − 05
202935_s_at
SOX9
SRY-box 9
5,971114
3,30E − 05
202023_at
EFNA1
Ephrin-A1
6,164569
3,30E − 05
202393_s_at
KLF10
Kruppel-like factor 10
6,194552
0,000195
213146_at
KDM6B
Lysine demethylase 6B
6,203146
1,90E − 05
205193_at
MAFF
v-maf avian musculoaponeurotic fibrosarcoma oncogene homolog F
6,2941
2,00E − 06
209457_at
DUSP5
Dual specificity phosphatase 5
6,639157
1,30E − 05
206029_at
ANKRD1
Ankyrin repeat domain 1
6,65759
0,008591
209283_at
CRYAB
Crystallin alpha B
6,703897
0,000118
201693_s_at
EGR1
Early growth response 1
7,056731
4,10E − 05
212099_at
RHOB
ras homolog family member B
7,300524
0,000406
219682_s_at
TBX3
T-box 3
7,722136
5,80E − 05
201473_at
JUNB
jun B proto-oncogene
8,322402
7,00E − 06
200664_s_at
DNAJB1
DnaJ heat shock protein family (Hsp40) member B1
8,586082
2,00E − 05
205828_at
MMP3
Matrix metallopeptidase 3
8,711976
1,90E − 05
201169_s_at
BHLHE40
Basic helix-loop-helix family member e40
8,870405
0,00011
203665_at
HMOX1
Heme oxygenase 1
9,32433
0,000544
202643_s_at
TNFAIP3
TNF alpha induced protein 3
9,573192
2,50E − 05
205207_at
IL6
Interleukin 6
10,18236
3,00E − 06
202388_at
RGS2
Regulator of G-protein signaling 2
10,25318
1,40E − 05
204472_at
GEM
GTP binding protein overexpressed in skeletal muscle
10,8003
1,00E − 06
202149_at
NEDD9
Neural precursor cell expressed, developmentally down-regulated 9
11,06553
2,50E − 05
219480_at
SNAI1
Snail family zinc finger 1
11,70457
2,00E − 06
218839_at
HEY1
hes related family bHLH transcription factor with YRPW motif 1
12,07541
6,00E − 06
206115_at
EGR3
Early growth response 3
14,19194
1,20E − 05
204470_at
CXCL1
C-X-C motif chemokine ligand 1
17,61827
2,00E − 06
204621_s_at
NR4A2
Nuclear receptor subfamily 4 group A member 2
18,77837
0
209774_x_at
CXCL2
C-X-C motif chemokine ligand 2
19,02731
9,00E − 06
202859_x_at
CXCL8
C-X-C motif chemokine ligand 8
19,89039
1,00E − 06
202340_x_at
NR4A1
Nuclear receptor subfamily 4 group A member 1
20,74943
1,00E − 06
209189_at
FOS
FBJ murine osteosarcoma viral oncogene homolog
23,36051
1,00E − 06
202672_s_at
ATF3
Activating transcription factor 3
24,18432
0
202768_at
FOSB
FBJ murine osteosarcoma viral oncogene homolog B
32,92245
0
207978_s_at
NR4A3
Nuclear receptor subfamily 4 group A member 3
43,80428
1,00E − 06
117_at
HSPA6
Heat shock protein family A (Hsp70) member 6
90,82389
0Figure 2
Transcriptional profiles of BEAS-2B cells infected with RSV for 4 h. (a) Hierarchical clustering of differentially expressed genes from dataset GSE3397. (b) Hierarchical clustering of differentially expressed genes from dataset GSE6802. RowZ-scores were calculated based on normalized expression data. The colors from green to red represent the transition of decreased to increased expression.
(a)
(b)


### 3.3. Functional Analysis

To obtain a biological interpretation of the transcriptional signature of RSV-infected BEAS-2B cells and compare with those reported by previous studies, enrichment analysis was performed with the online platform for functional analysis InnateDB [21]. Based on a FDR < 0.1, DEGs identified for GSE3397 were enriched in pathways related to Chromatin organization, histone acetylation, signaling by NOTCH, IL1, Integrin-linked kinase signaling, EPO signaling pathway, VEGF signaling pathway, platelet degranulation, p73 transcription factor network, IL-7 signaling, p53 signaling pathway, and others (Figure 3(a) and Supplementary Data 1). Of interest, Huang and collaborators (2008) [13] reported gene overrepresentation within p53 signaling pathway, but only after 24 h following RSV infection of BEAS-2B cells. After 4 h following RSV infection, Huang and collaborators (2008) [13] only found a significant association with neuroactive ligand-receptor interaction pathway, which was not overrepresented in the present analysis. In contrast, DEGs resultant from dataset GSE6802 were enriched in pathways related to AP-1 transcription factor, ATF-2 transcription factor, IL-6 signaling, SMAD function, signaling by TGFBR, HIF-1α transcription factor, signaling by CD40/CD40L, signaling by MAPK, signaling by innate immune receptors, and others (Figure 3(b) and Supplementary Data 1). Some of those pathways as CD40 signaling are indeed commonly induced by a variety of viral respiratory infections [22], whereas several of those pathways could indicate novel directions for studying the host response against RSV. Six pathways were enriched by DEGs from both datasets, the EPO signaling pathway, FBXW7 Mutants and NOTCH1 in Cancer, IL1, p53 signaling pathway, p73 transcription factor network, and signaling by NOTCH1. The erythropoietin (EPO) gene is a primary target of HIF-1α transcription factor, whereas binding of HIF-1α to the EPO enhancer promoter region induces transcriptional programs that influence inflammation and infection processes [23]. In addition, expression of Dll4, a major NOTCH ligand, is upregulated in dendritic cells infected with RSV, whereas blockage of Dll4in vivo increased hyperreactivity of airways and mucus secretion that impacted the pathology of the disease, showing a key role of signaling by NOTCH in the regulation of immunity against RSV [24]. Moreover, besides modulations of the p53 signaling pathway by infection of RSVin vitro [10, 13], this pathway was found to be upregulated in whole blood of children with lower respiratory tract infection by RSV [25]. Taken together, those data point to key pathways which can impact infections of human bronchial epithelial cells with RSV.Figure 3
Pathway enrichment analysis with InnateDB. Differentially expressed genes from (a) GSE3397 or (b) GSE6802 were evaluated for overrepresentation in pathways annotated in databases as INOH, KEGG, NETPATH, PID NIC, and REACTOME.
(a)
(b)


### 3.4. Meta-Analysis Based Biomarker Signature of RSV-Infected BEAS-2B Cells

To determine a unique transcriptional signature of BEAS-2B cells induced by early infection with RSV, common DEGs for both datasets were further identified. The analysis retrieved a list of seventeen common genes:ABCC4,ARMC8,BCLAF1,EZH1,FAM118A,FAM208B,FUS,HSPH1,KAZN,MAP3K2,N6AMT1,PRMT2,S100PBP,SERPINA1,TLK2,ZNF322, andZNF337 (Figure 4). Despite particular features in expression data from both datasets, unsupervised hierarchical clustering analysis based on this signature revealed the formation of robust clusters between RSV-infected or uninfected BEAS-2B cells (Figure 4). Of note, human airway epithelial cells were shown to express ABCC4/MRP4, a transporter for uric acid and cAMP [26]. Mucosal production of uric acid was recently linked to particulate matter-induced allergic sensitization [26]; therefore RSV infection could trigger such a response and contribute to the development and severity of allergic responses to particulate matter [27]. Moreover, both ABCC4 and SERPINA1 are annotated into the platelet degranulation pathway (Figure 3(a)), suggesting a role in antiviral mechanisms from bronchial epithelial cells. After an initial encounter with RSV, the transcriptional activity of human bronchial epithelial cells is reprogrammed to counteract viruses and other pathogens [10], whereasMAP3K2 andZNF322 are clearly involved on the activation and regulation of MAP kinase signaling pathway [28, 29]. Indeed, RSV infection leads to the activation of p38 MAPK [30] and c-JUN kinase pathway, which negatively regulates the production of TNF-α in human epithelial cells [31] and might contribute to virus evasion from an early immune response. Interestingly, the biosignature also included BCLAF1, a molecule involved in processes as apoptosis, transcription and processing of RNA, and export of mRNA from the nucleus [32]. However, this nuclear protein was also implicated as a viral restriction factor targeted to degradation by human cytomegalovirus [32]. Moreover, EZH1 was shown to be involved in the methylation of histone 3 at lysine 27 (H3K27) of the HIV provirus in resting cells [33] and could thus exert a significant function in infections with RSV, whereby other genes such as N6AMT1, FUS, and PRMT2 are also involved in protein methylation. Indeed, using coimmunoprecipitation and mass spectrometry, recent work demonstrated that RSV nucleoprotein (N) interacts with protein arginine N-methyltransferase 5 (PRMT5) [34], suggesting that PRMT2 could also interact with RSV proteins and play an important role during infections of human bronchial epithelial cells. Several of the genes identified in this study have been poorly studied in the context of RSV infection, whereby none of them was previously reported as a biomarker of infections by this virus. Of note, except forFAM208B andKAZN, analysis conducted by Smith and collaborators (2012) [22] which included both datasets (GSE3397 and GSE6802) also identified the significant modulation of the genes included in the biomarker signature identified herein.Figure 4
Biomarker signature of BEAS-2B cells infected with RSV for 4 h. Hierarchical clustering of expression data forABCC4,ARMC8,BCLAF1,EZH1,FAM118A,FAM208B,FUS,HSPH1,KAZN,MAP3K2,N6AMT1,PRMT2,S100PBP,SERPINA1,TLK2,ZNF322, andZNF337 from (a) dataset GSE3397 and (b) dataset GSE6802. Row Z-scores were calculated based on normalized expression data. The colors from blue to red represent the transition of decreased to increased expression.
(a)
(b)


## 3.1. Dataset Selection and Preprocessing Analysis

To define a robust transcriptional signature of BEAS-2B acutely infected with RSV, two publicly available datasets, GSE3397 and GSE6802, were used to conduct a meta-analysis from which data were extracted for BEAS-2B cells infected with RSV for 4 h and controls. First, background subtracted expression data from GSE3397 (Figure1(a)) were preprocessed and normalized (Figure 1(b)). However, in a first attempt to conduct differential gene expression analysis usinglimma [19], there were no statistically significant differences in gene expression. Therefore, principal component analysis (PCA) was used to evaluate the expression profiles of each array and, except for arrays named here Control2 and RSV2, the consistent pattern of clustering in Figure 1(c) suggests a batch effect. After normalization, this effect was even more evident (Figure 1(d)), which led to the speculation that Huang and collaborators (2008) [13] analyzed only three microarray experiments from this dataset based on the assumption that differences found for those microarrays were due to failures in experimental procedures; however they did not consider or correct for batch effects. In view of those facts, the datasets were adjusted with Combat function for R, which removed such effects from GSE3397 expression data (Figure 1(e)). Batch correction of GSE3397 did not change the profiles of arrays Control2 and RSV2; nevertheless, those arrays were included in further analysis because the variation observed in this experiment could have a substantial impact over the final result. Even adverse experimental variations that may change the overall expression patterns of a dataset could be useful to power up the identification of genes that are robustly modulated in BEAS-2B cells infected with RSV. The expression dataset GSE6802 (Figure 1(f)) was also included in the analysis. PCA from expression data extracted from GEO demonstrates that most of the variability between the arrays is explained (76.6%) by the infection with RSV, as the standardized PC1 separates RSV-infected from control arrays (Figure 1(g)), whereas standardized PC2 (11.4%) separates one pair of arrays (RSV_3 and ctrl2) and, although these arrays are supposedly from different batches, clustering features of this axis also suggested a batch effect (Figure 1(g)). log2 transformation of data impacted the profile of array RSV_1 however did not change the profiles from RSV_3 and ctrl_2 (Figure 1(h)). Combat( ) function was also applied to the expression dataset GSE6802; however, PCA shows that the adjustment did not to improve further clustering between specific arrays (Supplementary Figure  1; see Supplementary Material available online at http://dx.doi.org/10.1155/2016/3605302). In view of that, downstream analyses were carried out with normalized log2 transformed data.Figure 1
Preprocessing analysis of GEO datasets GSE3397 and GSE6802. (a) Boxplot of GSE3397, raw expression data. (b) Boxplot of GSE3397, normalized expression data. (c) Principal component analysis of GSE3397, raw expression data. (d) Principal component analysis of GSE3397, normalized expression data. (e) Principal component analysis of GSE3397, normalized and batch corrected expression data. (f) Boxplot of GSE6802, RMA normalized expression data. (g) Principal component analysis of GSE6802, normalized expression data. (h) Principal component analysis of GSE6802, log2 transformed RMA normalized expression data.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 3.2. Differential Gene Expression

Next, linear model-based statistical analyses with a FDR < 0.05 were conducted to identify differentially expressed genes (DEGs). The dataset GSE3397 exhibited ninety-four DEGs (Figure2(a) and Table 1). Those genes are highly discordant from DEGs previously reported by Huang and collaborators (2008) [13], which identified 277 DEGs based on different statistical analysis and assumptions. Fifty genes were downregulated and forty-four were upregulated (Table 1). The differences found in this study might reflect the inclusion of all microarray experiments from controls and 4 h after RSV infection; exclusion of expression data from 24 h after RSV infection; distinct preprocessing approaches as normalizing method and batch effect correction; and the assessment of statistical significance with a linear model-based method and corrected P values. In contrast, 1965 DEGs were identified for the dataset GSE6802. The top hundred DEGs ranked by fold changes (Figure 2(b) and Table 2) included genes such asJUNB,KLF4,CXCL1,CXCL2, andIL6, which are in agreement with those reported by Mayer and collaborators (2007) [10]. Several factors should account for the notable differences in expression analysis from both datasets. First, different RSV strains were used to stimulate BEAS-2B cells. Second, experimental conditions of controls were also different, as control experiments from GSE3397 were incubated with vehicle (not specified) and those from GSE6802 were not stimulated. Third, despite both datasets being generated with affymetrix microarray platform, those include distinct versions, HU133 plus 2.0 for GSE3397 and HU133A 2.0 for GSE6802.Table 1
Differentially expressed genes identified in dataset GSE3397.
ProbeID
Gene symbol
Gene name
l
o
g
2 fold change
FDR
1560754_at
CMTM7
CKLF like MARVEL transmembrane domain containing 7
−1,54756
0,017104
239439_at
AFF4
AF4/FMR2 family member 4
−1,53581
0,023832
238929_at
SRSF8
Serine/arginine-rich splicing factor 8
−1,51887
0,018433
223142_s_at
UCK1
Uridine-cytidine kinase 1
−1,47939
0,017104
242636_at
PRCP
Prolylcarboxypeptidase
−1,45095
0,034358
228007_at
CEP85L
Centrosomal protein 85 kDa-like
−1,4103
0,017104
235573_at
HSPH1
Heat shock protein family H (Hsp110) member 1
−1,39959
0,0371
228391_at
CYP4V2
Cytochrome P450 family 4 subfamily V member 2
−1,38799
0,01671
219376_at
ZNF322
Zinc finger protein 322
−1,3491
0,046761
1553689_s_at
METTL6
Methyltransferase like 6
−1,34723
0,017104
242837_at
SRSF4
Serine/arginine-rich splicing factor 4
−1,34071
0,044693
237215_s_at
TFRC
Transferrin receptor
−1,32685
0,017104
208819_at
RAB8A
RAB8A, member RAS oncogene family
−1,32593
0,042264
236665_at
CCDC18
Coiled-coil domain containing 18
−1,31494
0,034201
206147_x_at
SCML2
Sex comb on midleg-like 2 (Drosophila)
−1,30586
0,016454
229325_at
ZZZ3
Zinc finger ZZ-type containing 3
−1,30495
0,017104
1565716_at
FUS
FUS RNA binding protein
−1,29415
0,049505
205062_x_at
ARID4A
AT-rich interaction domain 4A
−1,28877
0,033039
1552312_a_at
MFAP3
Microfibrillar associated protein 3
−1,28521
0,046511
223223_at
ARV1
ARV1 homolog, fatty acid homeostasis modulator
−1,27987
0,023832
232001_at
PRKCQ-AS1
PRKCQ antisense RNA 1
−1,27987
0,035983
233195_at
DNAI1
Dynein axonemal intermediate chain 1
−1,25963
0,047083
219094_at
ARMC8
Armadillo repeat containing 8
−1,25527
0,043392
235232_at
GMEB1
Glucocorticoid modulatory element binding protein 1
−1,2492
0,046511
218643_s_at
CRIPT
CXXC repeat containing interactor of PDZ3 domain
−1,24229
0,0371
1566851_at
TRIM42
Tripartite motif containing 42
−1,24057
0,042149
221821_s_at
KANSL2
KAT8 regulatory NSL complex subunit 2
−1,23799
0,017104
244115_at
FAM126A
Family with sequence similarity 126 member A
−1,23114
0,033039
215541_s_at
DIAPH1
Diaphanous related formin 1
−1,22774
0,033039
203196_at
ABCC4
ATP binding cassette subfamily C member 4
−1,22519
0,033039
225024_at
RPRD1B
Regulation of nuclear pre-mRNA domain containing 1B
−1,22264
0,043765
37860_at
ZNF337
Zinc finger protein 337
−1,22095
0,023832
212997_s_at
TLK2
Tousled like kinase 2
−1,21841
0,04814
225690_at
CDK12
Cyclin-dependent kinase 12
−1,21083
0,0371
232103_at
BPNT1
3′(2′), 5′-Bisphosphate nucleotidase 1
−1,20748
0,0371
224848_at
CDK6
Cyclin-dependent kinase 6
−1,20247
0,0371
214962_s_at
NUP160
Nucleoporin 160 kDa
−1,20247
0,046319
219629_at
FAM118A
Family with sequence similarity 118 member A
−1,19831
0,028374
212290_at
SLC7A1
Solute carrier family 7 member 1
−1,19748
0,042264
227187_at
CBLL1
Cbl proto-oncogene like 1, E3 ubiquitin protein ligase
−1,19582
0,030047
233208_x_at
CPSF2
Cleavage and polyadenylation specific factor 2
−1,19334
0,046319
230566_at
MORC2-AS1
MORC2 antisense RNA 1
−1,17691
0,0371
238795_at
FAM208B
Family with sequence similarity 208 member B
−1,17609
0,0371
204980_at
CLOCK
Clock circadian regulator
−1,17283
0,0371
238653_at
LRIG2
Leucine-rich repeats and immunoglobulin like domains 2
−1,17202
0,048527
229939_at
ENDOV
Endonuclease V
−1,16878
0,041349
218185_s_at
ARMC1
Armadillo repeat containing 1
−1,16151
0,046319
201083_s_at
BCLAF1
BCL2 associated transcription factor 1
−1,15509
0,049505
227840_at
C2orf76
Chromosome 2 open reading frame 76
−1,15109
0,042264
201686_x_at
API5
Apoptosis inhibitor 5
−1,14076
0,046761
221699_s_at
DDX50
DEAD-box helicase 50
1,140764
0,046511
1556178_x_at
TAF8
TATA-box binding protein associated factor 8
1,159096
0,034358
205623_at
ALDH3A1
Aldehyde dehydrogenase 3 family member A1
1,163927
0,049505
212495_at
KDM4B
Lysine demethylase 4B
1,193336
0,044693
1569057_s_at
MIA3
Melanoma inhibitory activity family member 3
1,193336
0,047866
222494_at
FOXN3
Forkhead box N3
1,19582
0,048527
223311_s_at
MTA3
Metastasis associated 1 family member 3
1,19582
0,041439
215424_s_at
SNW1
SNW domain containing 1
1,196649
0,049505
213478_at
KAZN
Kazrin, periplakin interacting protein
1,19914
0,025143
227864_s_at
MVB12A
Multivesicular body subunit 12A
1,201636
0,030287
228674_s_at
EML4
Echinoderm microtubule associated protein like 4
1,204137
0,040345
224196_x_at
DPH5
Diphthamide biosynthesis 5
1,205808
0,025143
224652_at
CCNY
Cyclin Y
1,207481
0,046761
212968_at
RFNG
RFNG O-fucosylpeptide 3-beta-N-acetylglucosaminyltransferase
1,211673
0,0371
1555486_a_at
PRR5L
Proline rich 5 like
1,212513
0,017104
232837_at
KIF13A
Kinesin family member 13A
1,214195
0,042264
224320_s_at
MCM8
Minichromosome maintenance 8 homologous recombination repair factor
1,217566
0,033039
230131_x_at
ARSD
Arylsulfatase D
1,221793
0,0371
218225_at
ECSIT
ECSIT signalling integrator
1,224336
0,034358
222610_s_at
S100PBP
S100P binding protein
1,226885
0,030047
32259_at
EZH1
Enhancer of zeste 1 polycomb repressive complex 2 subunit
1,229439
0,0371
203854_at
CFI
Complement factor I
1,232852
0,042264
221600_s_at
AAMDC
Adipogenesis associated, Mth938 domain containing
1,260503
0,0371
209558_s_at
HIP1R
Huntingtin interacting protein 1 related
1,263127
0,042264
224814_at
DPP7
Dipeptidyl peptidase 7
1,26488
0,016454
232280_at
SLC25A29
Solute carrier family 25 member 29
1,277214
0,030047
228424_at
NAALADL1
N-Acetylated alpha-linked acidic dipeptidase-like 1
1,286989
0,042264
203409_at
DDB2
Damage specific DNA binding protein 2
1,288775
0,023832
229975_at
BMPR1B
Bone morphogenetic protein receptor type 1B
1,297739
0,034358
227073_at
MAP3K2
Mitogen-activated protein kinase kinase kinase 2
1,297739
0,017104
225347_at
ARL8A
ADP ribosylation factor like GTPase 8A
1,298639
0,02672
221774_x_at
SUPT20H
SPT20 homolog, SAGA complex component
1,308578
0,016454
223679_at
CTNNB1
Catenin beta 1
1,318594
0,018043
227679_at
HDAC11
Histone deacetylase 11
1,328686
0,044693
220020_at
XPNPEP3
X-Prolyl aminopeptidase 3, mitochondrial
1,342573
0,031097
203199_s_at
MTRR
5-Methyltetrahydrofolate-homocysteine methyltransferase reductase
1,360371
0,017104
228722_at
PRMT2
Protein arginine methyltransferase 2
1,370783
0,016454
228951_at
SLC38A7
Solute carrier family 38 member 7
1,431969
0,016454
217529_at
ORAI2
ORAI calcium release-activated calcium modulator 2
1,453973
0,043775
220311_at
N6AMT1
N-6 adenine-specific DNA methyltransferase 1 (putative)
1,460032
0,017104
213402_at
ZNF787
Zinc finger protein 787
1,469169
0,017104
226055_at
ARRDC2
Arrestin domain containing 2
1,477338
0,017104
219756_s_at
POF1B
Premature ovarian failure, 1B
1,580083
0,016454
202833_s_at
SERPINA1
Serpin peptidase inhibitor, clade A (alpha-1 antiproteinase, antitrypsin), and member 1
2,488023
0,0371Table 2
Top hundred differentially expressed genes identified in dataset GSE6802.
ProbeID
Gene symbol
Gene name
l
o
g
2 fold change
FDR
212615_at
CHD9
Chromodomain helicase DNA binding protein 9
−3,69609
0,00131
221840_at
PTPRE
Protein tyrosine phosphatase, receptor type E
−3,56524
0,000195
220817_at
TRPC4
Transient receptor potential cation channel subfamily C member 4
−3,39168
0,001582
221703_at
BRIP1
BRCA1 interacting protein C-terminal helicase 1
−2,88786
0,021463
207012_at
MMP16
Matrix metallopeptidase 16
−2,82647
0,000119
219494_at
RAD54B
RAD54 homolog B (S. cerevisiae)
−2,81279
0,000177
207034_s_at
GLI2
GLI family zinc finger 2
−2,79723
0,005157
203518_at
LYST
Lysosomal trafficking regulator
−2,75872
5,90E − 05
205282_at
LRP8
LDL receptor related protein 8
−2,7549
0,000311
214440_at
NAT1
N-Acetyltransferase 1 (arylamine N-acetyltransferase)
−2,68515
0,001777
219627_at
ZNF767P
Zinc finger family member 767, pseudogene
−2,67957
0,00024
218984_at
PUS7
Pseudouridylate synthase 7 (putative)
−2,67586
0,001308
206554_x_at
SETMAR
SET domain and mariner transposase fusion gene
−2,63536
0,002432
219779_at
ZFHX4
Zinc finger homeobox 4
−2,62624
0,001411
213103_at
STARD13
StAR related lipid transfer domain containing 13
−2,57219
0,002525
210138_at
RGS20
Regulator of G-protein signaling 20
−2,55974
0,000415
204291_at
ZNF518A
Zinc finger protein 518A
−2,54383
9,70E − 05
204651_at
NRF1
Nuclear respiratory factor 1
−2,49147
0,003659
205408_at
MLLT10
Myeloid/lymphoid or mixed-lineage leukemia; translocated to, 10
−2,48975
5,10E − 05
219581_at
TSEN2
tRNA splicing endonuclease subunit 2
−2,45377
0,001774
218242_s_at
SUV420H1
Lysine methyltransferase 5B
−2,44698
0,000754
203242_s_at
PDLIM5
PDZ and LIM domain 5
−2,43851
0,001699
203868_s_at
VCAM1
Vascular cell adhesion molecule 1
−2,43513
0,000761
220206_at
ZMYM1
Zinc finger MYM-type containing 1
−2,36362
0,008439
207616_s_at
TANK
TRAF family member associated NFKB activator
−2,34567
0,000424
218303_x_at
KRCC1
Lysine-rich coiled-coil 1
−2,34567
0,003187
218490_s_at
ZNF302
Zinc finger protein 302
−2,32785
0,001816
206876_at
SIM1
Single-minded family bHLH transcription factor 1
−2,32624
0,001681
219128_at
C2orf42
Chromosome 2 open reading frame 42
−2,28628
0,002926
212861_at
MFSD5
Major facilitator superfamily domain containing 5
−2,27048
0,000823
218653_at
SLC25A15
Solute carrier family 25 member 15
−2,25636
0,000562
206943_at
TGFBR1
Transforming growth factor beta receptor I
−2,24856
0,025349
201995_at
EXT1
Exostosin glycosyltransferase 1
−2,247
0,000421
221430_s_at
RNF146
Ring finger protein 146
−2,23457
0,001084
212286_at
ANKRD12
Ankyrin repeat domain 12
−2,2253
0,00029
219544_at
BORA
Bora, aurora kinase A activator
−2,21914
0,000333
210455_at
R3HCC1L
R3H domain and coiled-coil containing 1 like
−2,2176
0,0039
219459_at
POLR3B
Polymerase (RNA) III subunit B
−2,2176
0,000832
219078_at
GPATCH2
G-patch domain containing 2
−2,19923
0,000723
204547_at
RAB40B
RAB40B, member RAS oncogene family
−2,17648
0,001741
209760_at
KIAA0922
KIAA0922
−2,17347
0,001048
218791_s_at
KATNBL1
Katanin regulatory subunit B1 like 1
−2,17347
0,001187
205173_x_at
CD58
CD58 molecule
−2,17196
0,00022
204352_at
TRAF5
TNF receptor associated factor 5
−2,16895
0,002659
212441_at
KIAA0232
KIAA0232
−2,16595
0,006084
204236_at
FLI1
Fli-1 proto-oncogene, ETS transcription factor
−2,15397
0,005141
203072_at
MYO1E
Myosin IE
−2,15248
0,000154
219904_at
ZSCAN5A
Zinc finger and SCAN domain containing 5A
−2,14801
0,00144
219133_at
OXSM
3-Oxoacyl-ACP synthase, mitochondrial
−2,12285
0,002424
205798_at
IL7R
Interleukin 7 receptor
−2,11257
0,00506
205476_at
CCL20
C-C motif chemokine ligand 20
4,613942
9,50E − 05
213497_at
ABTB2
Ankyrin repeat and BTB domain containing 2
4,623547
1,40E − 05
219179_at
DACT1
Dishevelled-binding antagonist of beta-catenin 1
4,642816
9,00E − 06
219228_at
ZNF331
Zinc finger protein 331
4,723971
6,00E − 06
213139_at
SNAI2
Snail family zinc finger 2
4,76673
1,40E − 05
218177_at
CHMP1B
Charged multivesicular body protein 1B
4,806544
1,00E − 05
203304_at
BAMBI
BMP and activin membrane-bound inhibitor
4,826576
3,00E − 06
201631_s_at
IER3
Immediate early response 3
4,833271
3,00E − 06
218559_s_at
MAFB
v-maf avian musculoaponeurotic fibrosarcoma oncogene homolog B
4,870264
0,000468
220266_s_at
KLF4
Kruppel-like factor 4 (gut)
4,890561
0,00022
209211_at
KLF5
Kruppel-like factor 5 (intestinal)
4,924578
0,002036
209681_at
SLC19A2
Solute carrier family 19 member 2
4,927992
5,90E − 05
205266_at
LIF
Leukemia inhibitory factor
4,955395
2,20E − 05
204790_at
SMAD7
SMAD family member 7
5,073566
0,000283
221667_s_at
HSPB8
Heat shock protein family B (small) member 8
5,422657
2,90E − 05
212665_at
TIPARP
TCDD-inducible poly(ADP-ribose) polymerase
5,525098
1,00E − 05
202935_s_at
SOX9
SRY-box 9
5,971114
3,30E − 05
202023_at
EFNA1
Ephrin-A1
6,164569
3,30E − 05
202393_s_at
KLF10
Kruppel-like factor 10
6,194552
0,000195
213146_at
KDM6B
Lysine demethylase 6B
6,203146
1,90E − 05
205193_at
MAFF
v-maf avian musculoaponeurotic fibrosarcoma oncogene homolog F
6,2941
2,00E − 06
209457_at
DUSP5
Dual specificity phosphatase 5
6,639157
1,30E − 05
206029_at
ANKRD1
Ankyrin repeat domain 1
6,65759
0,008591
209283_at
CRYAB
Crystallin alpha B
6,703897
0,000118
201693_s_at
EGR1
Early growth response 1
7,056731
4,10E − 05
212099_at
RHOB
ras homolog family member B
7,300524
0,000406
219682_s_at
TBX3
T-box 3
7,722136
5,80E − 05
201473_at
JUNB
jun B proto-oncogene
8,322402
7,00E − 06
200664_s_at
DNAJB1
DnaJ heat shock protein family (Hsp40) member B1
8,586082
2,00E − 05
205828_at
MMP3
Matrix metallopeptidase 3
8,711976
1,90E − 05
201169_s_at
BHLHE40
Basic helix-loop-helix family member e40
8,870405
0,00011
203665_at
HMOX1
Heme oxygenase 1
9,32433
0,000544
202643_s_at
TNFAIP3
TNF alpha induced protein 3
9,573192
2,50E − 05
205207_at
IL6
Interleukin 6
10,18236
3,00E − 06
202388_at
RGS2
Regulator of G-protein signaling 2
10,25318
1,40E − 05
204472_at
GEM
GTP binding protein overexpressed in skeletal muscle
10,8003
1,00E − 06
202149_at
NEDD9
Neural precursor cell expressed, developmentally down-regulated 9
11,06553
2,50E − 05
219480_at
SNAI1
Snail family zinc finger 1
11,70457
2,00E − 06
218839_at
HEY1
hes related family bHLH transcription factor with YRPW motif 1
12,07541
6,00E − 06
206115_at
EGR3
Early growth response 3
14,19194
1,20E − 05
204470_at
CXCL1
C-X-C motif chemokine ligand 1
17,61827
2,00E − 06
204621_s_at
NR4A2
Nuclear receptor subfamily 4 group A member 2
18,77837
0
209774_x_at
CXCL2
C-X-C motif chemokine ligand 2
19,02731
9,00E − 06
202859_x_at
CXCL8
C-X-C motif chemokine ligand 8
19,89039
1,00E − 06
202340_x_at
NR4A1
Nuclear receptor subfamily 4 group A member 1
20,74943
1,00E − 06
209189_at
FOS
FBJ murine osteosarcoma viral oncogene homolog
23,36051
1,00E − 06
202672_s_at
ATF3
Activating transcription factor 3
24,18432
0
202768_at
FOSB
FBJ murine osteosarcoma viral oncogene homolog B
32,92245
0
207978_s_at
NR4A3
Nuclear receptor subfamily 4 group A member 3
43,80428
1,00E − 06
117_at
HSPA6
Heat shock protein family A (Hsp70) member 6
90,82389
0Figure 2
Transcriptional profiles of BEAS-2B cells infected with RSV for 4 h. (a) Hierarchical clustering of differentially expressed genes from dataset GSE3397. (b) Hierarchical clustering of differentially expressed genes from dataset GSE6802. RowZ-scores were calculated based on normalized expression data. The colors from green to red represent the transition of decreased to increased expression.
(a)
(b)


## 3.3. Functional Analysis

To obtain a biological interpretation of the transcriptional signature of RSV-infected BEAS-2B cells and compare with those reported by previous studies, enrichment analysis was performed with the online platform for functional analysis InnateDB [21]. Based on a FDR < 0.1, DEGs identified for GSE3397 were enriched in pathways related to Chromatin organization, histone acetylation, signaling by NOTCH, IL1, Integrin-linked kinase signaling, EPO signaling pathway, VEGF signaling pathway, platelet degranulation, p73 transcription factor network, IL-7 signaling, p53 signaling pathway, and others (Figure 3(a) and Supplementary Data 1). Of interest, Huang and collaborators (2008) [13] reported gene overrepresentation within p53 signaling pathway, but only after 24 h following RSV infection of BEAS-2B cells. After 4 h following RSV infection, Huang and collaborators (2008) [13] only found a significant association with neuroactive ligand-receptor interaction pathway, which was not overrepresented in the present analysis. In contrast, DEGs resultant from dataset GSE6802 were enriched in pathways related to AP-1 transcription factor, ATF-2 transcription factor, IL-6 signaling, SMAD function, signaling by TGFBR, HIF-1α transcription factor, signaling by CD40/CD40L, signaling by MAPK, signaling by innate immune receptors, and others (Figure 3(b) and Supplementary Data 1). Some of those pathways as CD40 signaling are indeed commonly induced by a variety of viral respiratory infections [22], whereas several of those pathways could indicate novel directions for studying the host response against RSV. Six pathways were enriched by DEGs from both datasets, the EPO signaling pathway, FBXW7 Mutants and NOTCH1 in Cancer, IL1, p53 signaling pathway, p73 transcription factor network, and signaling by NOTCH1. The erythropoietin (EPO) gene is a primary target of HIF-1α transcription factor, whereas binding of HIF-1α to the EPO enhancer promoter region induces transcriptional programs that influence inflammation and infection processes [23]. In addition, expression of Dll4, a major NOTCH ligand, is upregulated in dendritic cells infected with RSV, whereas blockage of Dll4in vivo increased hyperreactivity of airways and mucus secretion that impacted the pathology of the disease, showing a key role of signaling by NOTCH in the regulation of immunity against RSV [24]. Moreover, besides modulations of the p53 signaling pathway by infection of RSVin vitro [10, 13], this pathway was found to be upregulated in whole blood of children with lower respiratory tract infection by RSV [25]. Taken together, those data point to key pathways which can impact infections of human bronchial epithelial cells with RSV.Figure 3
Pathway enrichment analysis with InnateDB. Differentially expressed genes from (a) GSE3397 or (b) GSE6802 were evaluated for overrepresentation in pathways annotated in databases as INOH, KEGG, NETPATH, PID NIC, and REACTOME.
(a)
(b)


## 3.4. Meta-Analysis Based Biomarker Signature of RSV-Infected BEAS-2B Cells

To determine a unique transcriptional signature of BEAS-2B cells induced by early infection with RSV, common DEGs for both datasets were further identified. The analysis retrieved a list of seventeen common genes:ABCC4,ARMC8,BCLAF1,EZH1,FAM118A,FAM208B,FUS,HSPH1,KAZN,MAP3K2,N6AMT1,PRMT2,S100PBP,SERPINA1,TLK2,ZNF322, andZNF337 (Figure 4). Despite particular features in expression data from both datasets, unsupervised hierarchical clustering analysis based on this signature revealed the formation of robust clusters between RSV-infected or uninfected BEAS-2B cells (Figure 4). Of note, human airway epithelial cells were shown to express ABCC4/MRP4, a transporter for uric acid and cAMP [26]. Mucosal production of uric acid was recently linked to particulate matter-induced allergic sensitization [26]; therefore RSV infection could trigger such a response and contribute to the development and severity of allergic responses to particulate matter [27]. Moreover, both ABCC4 and SERPINA1 are annotated into the platelet degranulation pathway (Figure 3(a)), suggesting a role in antiviral mechanisms from bronchial epithelial cells. After an initial encounter with RSV, the transcriptional activity of human bronchial epithelial cells is reprogrammed to counteract viruses and other pathogens [10], whereasMAP3K2 andZNF322 are clearly involved on the activation and regulation of MAP kinase signaling pathway [28, 29]. Indeed, RSV infection leads to the activation of p38 MAPK [30] and c-JUN kinase pathway, which negatively regulates the production of TNF-α in human epithelial cells [31] and might contribute to virus evasion from an early immune response. Interestingly, the biosignature also included BCLAF1, a molecule involved in processes as apoptosis, transcription and processing of RNA, and export of mRNA from the nucleus [32]. However, this nuclear protein was also implicated as a viral restriction factor targeted to degradation by human cytomegalovirus [32]. Moreover, EZH1 was shown to be involved in the methylation of histone 3 at lysine 27 (H3K27) of the HIV provirus in resting cells [33] and could thus exert a significant function in infections with RSV, whereby other genes such as N6AMT1, FUS, and PRMT2 are also involved in protein methylation. Indeed, using coimmunoprecipitation and mass spectrometry, recent work demonstrated that RSV nucleoprotein (N) interacts with protein arginine N-methyltransferase 5 (PRMT5) [34], suggesting that PRMT2 could also interact with RSV proteins and play an important role during infections of human bronchial epithelial cells. Several of the genes identified in this study have been poorly studied in the context of RSV infection, whereby none of them was previously reported as a biomarker of infections by this virus. Of note, except forFAM208B andKAZN, analysis conducted by Smith and collaborators (2012) [22] which included both datasets (GSE3397 and GSE6802) also identified the significant modulation of the genes included in the biomarker signature identified herein.Figure 4
Biomarker signature of BEAS-2B cells infected with RSV for 4 h. Hierarchical clustering of expression data forABCC4,ARMC8,BCLAF1,EZH1,FAM118A,FAM208B,FUS,HSPH1,KAZN,MAP3K2,N6AMT1,PRMT2,S100PBP,SERPINA1,TLK2,ZNF322, andZNF337 from (a) dataset GSE3397 and (b) dataset GSE6802. Row Z-scores were calculated based on normalized expression data. The colors from blue to red represent the transition of decreased to increased expression.
(a)
(b)


## 4. Conclusions

The combined analysis of distinct datasets from BEAS-2B cells infected with RSV retrieved intriguing results, whereby using powerful statistical methods and assumptions this study identified a new set of biomarkers of early infection with RSV composed by seventeen genes:ABCC4,ARMC8,BCLAF1,EZH1,FAM118A,FAM208B,FUS,HSPH1,KAZN,MAP3K2,N6AMT1,PRMT2,S100PBP,SERPINA1,TLK2,ZNF322, andZNF337. This transcriptional signature could be useful for the development of molecular diagnosis tools as well as future investigations of processes involved in host-pathogen interactions.



---

*Source: 3605302-2016-05-04.xml*