Distinct transcriptional modules in the peripheral blood mononuclear cells response to human respiratory syncytial virus or to human rhinovirus in hospitalized infants with bronchiolitis

Human respiratory syncytial virus (HRSV) is the main cause of bronchiolitis during the first year of life, when infections by other viruses, such as rhinovirus, also occur and are clinically indistinguishable from those caused by HRSV. In hospitalized infants with bronchiolitis, the analysis of gene expression profiles from peripheral blood mononuclear cells (PBMC) may be useful for the rapid identification of etiological factors, as well as for developing diagnostic tests, and elucidating pathogenic mechanisms triggered by different viral agents. In this study we conducted a comparative global gene expression analysis of PBMC obtained from two groups of infants with acute viral bronchiolitis who were infected by HRSV (HRSV group) or by HRV (HRV group). We employed a weighted gene co-expression network analysis (WGCNA) which allows the identification of transcriptional modules and their correlations with HRSV or HRV groups. This approach permitted the identification of distinct transcription modules for the HRSV and HRV groups. According to these data, the immune response to HRSV infection—comparatively to HRV infection—was more associated to the activation of the interferon gamma signaling pathways and less related to neutrophil activation mechanisms. Moreover, we also identified host-response molecular markers that could be used for etiopathogenic diagnosis. These results may contribute to the development of new tests for respiratory virus identification. The finding that distinct transcriptional profiles are associated to specific host responses to HRSV or to HRV may also contribute to the elucidation of the pathogenic mechanisms triggered by different respiratory viruses, paving the way for new therapeutic strategies.

Human respiratory syncytial virus (HRSV) is the main cause of bronchiolitis during the first year of life, when infections by other viruses, such as rhinovirus, also occur and are clinically indistinguishable from those caused by HRSV. In hospitalized infants with bronchiolitis, the analysis of gene expression profiles from peripheral blood mononuclear cells (PBMC) may be useful for the rapid identification of etiological factors, as well as for developing diagnostic tests, and elucidating pathogenic mechanisms triggered by different viral agents. In this study we conducted a comparative global gene expression analysis of PBMC obtained from two groups of infants with acute viral bronchiolitis who were infected by HRSV (HRSV group) or by HRV (HRV group). We employed a weighted gene co-expression network analysis (WGCNA) which allows the identification of transcriptional modules and their correlations with HRSV or HRV groups. This approach permitted the identification of distinct transcription modules for the HRSV and HRV groups. According to these data, the immune response to HRSV infection-comparatively to HRV infection-was more associated to the activation of the interferon gamma signaling pathways and less related to neutrophil activation mechanisms. Moreover, we also identified host-response molecular markers that could be used for etiopathogenic diagnosis. These results may contribute to the development of new tests for respiratory virus identification. The finding that distinct transcriptional profiles are associated to specific host responses to HRSV or to HRV may also contribute to the elucidation of the pathogenic mechanisms triggered by different respiratory viruses, paving the way for new therapeutic strategies. PLOS

Ethics statement
This study was approved by the Research Ethics Committee of the Hospital Universitário da Universidade de São Paulo under number 1011/10. The infants' legal guardians signed the Informed Consent Form after the presentation of the study by one of the authors and answered standardized questionnaires for the obtaining clinical, demographic and epidemiological data. All infants were submitted to molecular analysis of nasopharyngeal secretion for viral identification. Patients' peripheral blood samples were collected for genomic analysis.

Patient characteristics
A total of 12 out of 124 infants under 6 months of age-hospitalized between 2013 and 2015 at the Hospital Universitário da Universidade de São Paulo-were selected from a prospective cohort study on the etiology of acute viral bronchiolitis [15]. Diagnosis of bronchiolitis was defined as the first wheezing crisis, beginning no more than 3 days before hospital admission.
In this study we included only infants infected with HRSV A ON1 or HRV as a single agent. Patients' demographic and clinical data are listed in Table 1. The median age of the infants of the HRSV group and the HRV group were similar (41.5 days and 64 days, respectively; p = 0.8). There was a predominance of males in both groups (66.6% in the HRSV group and 83.3% in the HRV group). The mean length of hospital stay in the HRSV group was 6 days and in the HRV group was 2 days.

Sample collection and total RNA extraction
Whole blood was collected in an EDTA-containing tube within a maximum of 24 hours after hospital admission. The peripheral blood mononuclear cells (PBMC) were immediately separated by gradient centrifugation using Ficoll-Paque Plus (GE Heathcare, cat. no. 17144002). After cell separation, the PBMC were collected, preserved in RNAlater (Qiagen, cat. no. 76106)

Microarray hybridization and global gene expression
In order to determine gene expression profiles, 44 K DNA microarrays (Whole Human Genome Microarray Kit, Agilent Technologies, cat no. G4112F, Santa Clara, CA, USA) were used. The procedures for hybridization using the fluorescent dye Cy3 followed the manufacturer's protocols (One-Color Microarray-Based Gene Expression Analysis-Quick Amp Labeling). The microarray images were captured by the reader Agilent Bundle, according to the parameters recommended for bio-arrays and extracted by Agilent Feature Extraction software version 9.5.3. Spots with two or more flags (low intensity, saturation, controls, etc.) were considered as NA, that is, without valid expression value. The R software version 2.11.1 and an in-house script were used for: i) excluding transcript spots presenting one or more NAs; iii) converting gene expression values to log base 2. Through this procedure the gene expression matrix with only expressed transcripts were obtained. Data normalization was performed using R software and the Lowess method [16]. TMEV software version 4.6.1 and t-test with p<0.05 and fold-change value � 2.0 was used for obtaining the differentially expressed genes for the comparison (HRSV versus HRV). All microarray raw data have been deposited in Expression Omnibus (GEO) public database under accession number GSE124124.

Weighted gene coexpression network analysis (WGCNA)
WGCNA is a method that identifies and characterizes gene modules whose members share strong coexpression. Networks were constructed using the WGCNA R software package [17]. Pearson's correlation coefficient was used for obtaining gene coexpression similarity measures and for the subsequent construction of an adjacency matrix using soft power and topological overlap matrix (TOM). Soft-thresholding process transforms the correlation matrix to mimic the scale free topology. TOM is used to filter weak connections during network construction. Module identification is based on TOM and in average linkage hierarchical clustering. Finally, dynamic cut-tree algorithm was used for dendrogram's branch selection. The module eigengene (ME) is defined as the first principal component of a given module, which can be considered a representative of the gene expression profiles in a module. Module Membership (MM), also known as eigengene-based connectivity (kME), is defined as the correlation of each gene expression profile with the module eigengene of a given module [17]. Module-HRSV or HRV association. Firstly, we obtained the gene significance (GS), which is a value of the correlation between the trait (here is represented by HRSV or HRV groups) and the gene expression values. The mean GS for a particular module is considered as a measure of module significance (MS). The GS values were obtained using Pearson's correlation and to assign a p-value to the module significance, we used Student's t test. The modules, which presented high positive correlation value with HRSV or HRV groups (r � 0.6 and p < 0.05) were selected for biological functional analysis.
Intramodular analysis for hub selection. The MM and GS values were used for gene categorization. Genes presenting high GS and MM were considered as hubs in the module and significantly associated with HRSV or HRV groups. We plotted all gene values in a MM (xaxis) vs GS graphic (y-axis).

Functional enrichment analysis
KEGG pathway enrichment was performed using Enrichr online web-based tool [18] to analyze the network modules that are associated with the HRSV or HRV groups (p < 0.05). The biological function analysis for the selected genes was done using the Gene Ontology database.

Results
Transcriptional analysis revealed that 283 out of 6,615 GO annotated genes were differentially expressed (DE, fold-change � 2.0) when compared between the HRSV and HRV groups. In the HRSV group most of the DE genes (204 genes) were hyper-expressed and only 79 genes were hypo-expressed when compared with the HRV group. The gene CCDC177 (C14orf162) presented the highest fold-change value (16.3). Moreover, 22 genes were differentially expressed and presented high fold-change values between 16.3 and 4.0 (S1 Fig). Table 2 lists these high fold-change DE genes, where 21 genes were hyper-expressed and one was hypoexpressed in HRSV group.
The enrichment analysis of the DE genes, based on Biological Process (BP) terms of Gene Ontology database, showed that 99 DE genes (60 hyper-expressed and 39 hypo-expressed genes in the HRSV group in comparison with the HRV group) are significantly over- represented in BP terms. Twenty-one hyper-expressed genes are involved in immune response processes, such as inflammatory response, apoptosis, response to interferon-γ, and neutrophil chemotaxis. Thirteen hypo-expressed genes are also related to immune response, involving neutrophil degranulation, inflammatory response, apoptosis and phagocytosis (Fig 1; S1 and S2 Tables).

WGCNA
The coexpression network was constructed by Pearson's correlation between all genes and considering a soft-thresholding power β of 22, resulting in a scale-free network. Following dynamic tree cut, the hierarchical clustering dendrogram identified 17 distinct gene modules (Fig 2). Modules are defined as branches of the network dendrogram and identified by different color. Modules size ranged from 51 (grey60 module) to 1,196 (turquoise module) genes.
Genes not classified in any correlated module were grouped in a grey module. After the modules were generated, HRSV and HRV groups were then correlated with the modules. A total of five modules were significantly (p<0.05) associated with HRSV. Three of those five modules were positively associated with HRSV group: midnight blue (MS = 0.79, p = 0.002), turquoise (MS = 0.69, p = 0.01) and tan (MS = 0.63, p = 0.03) modules. Others two modules were positively associated with HRV group: purple (MS = 0.78, p = 0.003) and brown (MS = 0.59, p = 0.04) modules. None module was associated with gender ( Fig 2B).
The functional profiles of these five gene modules-based on enrichment analysis using KEGG pathways database-showed that the midnight blue and tan modules (positively associated to HRSV) contain a high proportion of genes (about 60%) involved in cell signaling and in immune response to viral infection. Conversely, the turquoise (positively associated to HRSV), brown, and purple modules (positively associated to HRV) have relatively few genes (26% in brown, 19% in turquoise, and 3.7% in purple) involved in cell/ immune response to viruses (Fig 3). A complete list of KEGG pathways found for these modules are listed in S3 Table. Then, we categorized the genes of the five modules positively associated with HRSV (midnight blue, turquoise, and tan) or with HRV (brown and purple) considering module membership (MM) and gene significant (GS) values for HRSV or HRV groups (S2 Fig). We identified a total of 42 and 57 genes with highest GS and MM values (GS and MM values � 0.75; p < 0.01 was considered significant; S4 and S5 Tables) for HRSV and for HRV groups, respectively. These genes were named here as HGS-hubs. The biological function of the HGS-hubs in HRSV and HRV groups is listed in Tables 3 and 4 In the HRSV group, 16 out of 46 HGS-hubs were related to immune and inflammatory responses, phagocytosis and apoptosis (Table 3). Moreover, the majority of these HGS-hubs is hyper-expressed in the HRSV group. It's interesting to note that 11 out of 28 HGS-hubs had fold-change values > 2.0. Seven of these genes are involved in immune responses: six are hyper-expressed (TPX2, CXCR6, DCST2, LGALS3BP, APOL1, and IDO1); and one (EFNB1) is hypo-expressed.
In the HRV group, comparatively to the HRSV group, just a small number (19 out of 71 genes) of HGS-hubs are involved in cell response to virus and majority of them is hyperexpressed in the HRV group (Table 4). Three hyper-expressed HGS-hubs are involved in immune response signaling pathways, such as negative regulation of NF-KappaB kinase signaling pathway (RASSF2 and CHIP) and TOR signaling pathway (CCDC88A). Only two HGShubs, CBL and PACS2, had fold-change values > 2.0 and are involved in inflammatory response and autophagic process, respectively.

Discussion
Bronchiolitis is frequently caused by HRSV but HRV is also an important etiological agent. Infections by these two viruses present similar clinical features, thus rendering etiological diagnosis difficult. Moreover, no vaccines or effective/specific therapies are available. The . Only the HRSV group is shown. Numbers inside each colored box are the correlation coefficients between the ME and the specific trait, with p-value between brackets. The same values are true for the HRV group, but with an opposite correlation coefficient signal. The more intense the color of the box, the more identification of the etiology in infant bronchiolitis by obtaining peripheral blood samples and performing molecular marker analyses would be of great utility in clinical practice. These methods are non-invasive and reflect quite well the host's immunological response [21], thus allowing a better understanding of the molecular mechanisms involved in the pathogenesis of HRSV infection.
Previous epidemiological and genotypical studies have indicated that the diseases caused by HRSV and HRV display some differences [4]. Infants, especially under 6 months of age, are more susceptible to hospitalization and death in HRSV infections [1,2]. Rodriguez-Fernandez et al [13] showed that infants with bronchiolitis presented different transcriptome profiles for distinct HRSV genotypes. For instance, the GA5 genotype was associated with diminished interferon activation and increased expression of genes involved in neutrophil activation, comparatively to other HRSV genotypes. Considering these facts, we selected infants under 6 months of age, who had a first episode of wheezing, and only infants with single infection with HRSV A ON1 genotype.
In this study, we found two distinct PBMC transcriptomic profiles and the differentially expressed genes subsets are related to different innate immune response of the patients with HRSV or HRV bronchiolitis. Mejias et al. [7] also showed a specific gene expression profile in infants infected with HRSV, characterized by overexpression of innate immunity and suppression of adaptive immunity, and associated with neutrophil response, inflammation, IFN activation and suppression of T and B cell genes. Here is important to recall that HRSV infects PBMCs during the immune response to viral challenge, when these cells are recruited to the respiratory trait, and differential immune response to HRSV is determined early after exposure [22].
Additionally, the leukogram profile (Table 1) showed significantly (Mann-Whitney test, p = 0.03) more leukocytes in infants with HRV infection than in infants with HRSV infection. Choi et al [23] also found leukocytosis in infants with HRV infection compared with other respiratory viruses (HRSV, influenza virus, parainfluenza virus, human coronavirus, human bocavirus, human metapneumovirus). However, all infants studied here had neutropenia (neutrophil/lymphocyte ratio � 1.0) and no significant difference between HRSV and HRV groups was found. It's interesting to mention that Choi et al [23] found a different neutrophil count, which neutrophil/lymphocyte ratio was over than 3.0 for HRV and less than 1.0 for HRSV infection. This neutropenia observed in both group of infants in our study probably was caused by neutrophil migration to lung as previously observed by Everard et al [24] analyzing bronchial secretion of children with HRSV infection.
It is well established that systemic neutrophil response correlates with increased disease severity and it is mediated by the neutrophil chemoattractant. In the fatal cases of HRSV inflammatory neutrophils were found widespread in lung tissue [25]. Moreover, Choi et al [23] showed that HRV is associated with neutrophil activation and inflammatory response while HRSV is associated with lymphocyte-mediated inflammation. In this study, the PBMC transcriptomic analysis revealed different expression profiles between HRSV and HRV infection. Five hypo-expressed genes in HRSV group-SYNGR1, AGL, LAMP2, ORM2, and SIRPB1 -are involved in neutrophil activation processes, while three genes expressing neutrophil chemoattractants-CCL15, CCL24, and LGALS3-were hyper-expressed in HRSV infection (S1 and S2 Tables). Our results indicated that the host response to HRSV, compared with HRV, are related to diminished neutrophil activation and increase neutrophil chemotaxis. This negatively (green) or positively (red) correlated is the module with the trait. Five modules presented significant association (p<0.05) with the HRSV group and these modules are indicated by black borders.  Distinct PBMC transcriptional modules in the response to HRSV or HRV in bronchiolitis  Distinct PBMC transcriptional modules in the response to HRSV or HRV in bronchiolitis response could well be associated to disease severity and, consequently, to longer periods of hospitalization (Table 1). Another subset of hyper-expressed genes-OAS2, PTAFR, PML, MEFV, CCL15, and CCL24-in HRSV group is related to interferon-gamma pathway (S1 Table). Interferon induces antiviral effectors encoding genes and exerts a protective effect [25,26]. Moreover, the balance between pro-and anti-inflammatory T-cell cytokine responses may determine the clinical outcome of HRSV infection. High levels of IFN-gamma might reduce viral replication but increase immunopathology. However, decreased IFN-gamma response has been reported to be associated with severe HRSV infection [27]. The treatment with systemic steroids is controversial regarding the efficacy of systemic steroid use in severe bronchiolitis. Pinto et al [28] showed that infants with severe illness had higher plasma cortisol levels than infants with mild disease, and a decreased IFN-gamma production by PBMCs in severely affected infants. Nevertheless, these two molecular mechanisms-neutrophil chemotaxis and IFN-gamma/cortisol-seem to be important therapeutic targets for modulating HRSVderived immunopathologies.
The gene coexpression network analysis revealed distinct modules (tan, midnight-blue, turquoise, purple and brown) containing genes involved in immune response, activation of interferon pathway and apoptosis. Module analysis identified HGS-hubs associated with HRSV and HRV groups. In network, hubs are important to maintain the structure of the network or modules [29][30][31]. Additionally, high gene significant (HGS) value implies that the gene is significantly associated to phenotypic/genotypic features of the group [17], i.e., HRSV or HRV infections and the gene is differentially expressed between groups. Therefore, these HGS-hubs might be used to differentiate these two etiologies.
Several HGS-hubs in both groups (14 out of 42 in HRSV group and 14 out of 57 in HRV group, Tables 3 and 4) are involved in immune response but the molecular mechanism differs between HRSV and HRV groups: the HGS-hubs involved in inflammatory response, apoptosis, and autophagy are not the same in each group. It's is interesting to mention that nine HGS-hubs had fold-change > 2.0 and are involved in immune response, inflammatory, and apoptosis. Moreover, eight of these genes might be potential therapeutic/vaccine targets and/ or molecular markers, as commented below. Distinct PBMC transcriptional modules in the response to HRSV or HRV in bronchiolitis Six out of eight HGS-hubs were previously described as potential therapeutic targets: i) LGALS3BP encodes for galectin 3 binding protein, an antiviral protein that interferes with human immunodeficiency virus type 1 replication diminishing the virions infectivity [32]; ii) EFNB1 codifies an ephrin-B1; members of Efns family are expressed in thymocytes and T cells and they are capable to modulate T cell responses and survival contributing to the integrity of an immune response [33]; iii) APOL1, codifies for an apolipoprotein 1, which is a mediator protein acting on IFN-activated genes [34,35] and the proteins of the apolipoprotein family, such as APOL6, has an HRSV antiviral activity related to apoptosis [36]; iv) IDO1 encoded protein was described to control viral infection by modulating specific metabolic events [37]. IDO1 was hyper-expressed after IFN-mediated induction in response to Influenza A viruses and this enzyme is involved in cytokine overproduction and T cell suppression [38]; v) CBL, alias C-CBL, encodes for RING finger E3 ubiquitin ligase and it is involved in antiviral process by activation IFN-I pathway [39]; and vi) PACS2, that codifies for a sorting protein PACS-2, interacts with a protein Nef of the human immunodeficiency virus type 1 [40]. This interaction mediates the MHC-I down-regulation and consequently the virus escape from host immunity [40].
Another strategy to combat HRSV could be a vaccine therapy. We identified one HGS-hub that might be a potent vaccine target: CXCR6, that codifies for a chemokine expressed by natural killer T (NKT) cells and it's essential for NKT cells trafficking [41]. It was described that the activation of NKT cells by intranasal coadministration of α-galactosylceramide, in a nasal vaccine against influenza, can potently enhance protective immune responses through increasing NKT cell population in nasal mucosa [41].
Finally, two HGS-hubs could be potential molecular markers for HRSV infection.
LGALS3BP, which is also a potential therapeutic target [32], and the other is DCST2 that codifies for a dendritic cell-specific transmembrane protein and its intracellular localization is due in response to toll-like receptor ligation [42]. These genes presented high fold-changes and were hyper-expressed in PBMC from infants infected with HRSV. Additionally, 22 DE genes had higher fold-changes (� 4.0) and six of them (one is the HGS-hub LGALS3BP) are involved in immune response, and other two genes related to transcriptional regulation are HGS-hubs (Table 2). Finally, the gene CCDC177 codifying a coiled-coil domain protein presented the highest fold-change 16.3.

Conclusions
The PBMC transcriptional profiles of hospitalized infants with HRSV or HRV bronchiolitis are different and probably correlated with distinctive etiopathogenic mechanisms. Additionally, the infants with HRSV had leukopenia and they usually demand more hospitalization days. The immune response to HRSV infection, comparatively to the response to HRV infection, appears to be more associated to the activation of the IFNγ signaling pathways and less related to neutrophil activation. Moreover, we also identified potential host-response molecular markers that could be used for HRSV or HRV etiopathogenic diagnosis. These results may contribute to the development of future tests for respiratory virus identification. Additionally, a better understanding of PBMC specific host responses to HRSV or HRV-here disclosed by different gene expression profiles-may serve to elucidate the pathogenic mechanisms triggered by different viral agents and, therefore, contribute to the development of new therapeutic approaches. Scatterplots between MM (xaxis) and GS (y-axis) of the genes in the i) midnight blue, turquoise, and tan modules of the HRSV group; ii) brown and purple modules of the HRV group. The red or blue dot lines indicate, respectively, the cut off of GS or MM values significantly for HRSV or HRV groups (p < 0.01). (TIF) S1