Upper Airway Cell Transcriptomics Identify a Major New Immunological Phenotype with Strong Clinical Correlates in Young Children with Acute Wheezing

This article is featured in In This Issue, p.1647

Exacerbations of asthma and wheeze are mostly triggered by respiratory viral infections and are one of the most common reasons for a child to seek emergency care (1). Previous studies from our group have demonstrated that rhinovirus (RV) species C is the most frequent viral pathogen detected in children who present to the local emergency department (ED) with an asthma exacerbation (2). However, the molecular mechanisms that determine susceptibility to RV and expression of respiratory symptoms are not well understood. Previous investigations of airway epithelial cells infected with RV in vitro found that type I and III IFN responses were deficient in adults with asthma, leading to impaired viral control and exaggerated secondary responses (3, 4). However, this finding was replicated in some studies but not others (5). Human adult asthmatic volunteers experimentally infected with RV in vivo have exaggerated IL-25 and IL-33 responses, which drive Th2 inflammation (6, 7). Although these data provide a plausible mechanism to link RV infection with the pathogenesis of asthma, they are based on experimental infections, and more data are needed in the context of natural exacerbations (1, 8). Notably, studies of naturally occurring virus-induced exacerbations report increased rather than decreased IFN responses (9–11). In this study, we have used an unbiased, systems biology approach to elucidate the innate immune mechanisms that are operating in the upper airways during natural exacerbations of asthma or wheezing in children (9, 10). Our findings provide new insights into the role of gene networks, particularly IRF7, and their relationship to clinical phenotypes in this disease.

The participants were part of an ongoing study examining the mechanisms of acute viral respiratory infection in children (MAVRIC). Cases included children aged 0–16 y presenting to the ED of a tertiary children’s hospital (Princess Margaret Hospital, Perth, WA, Australia) with acute wheeze and the availability of nasal fluid and swab specimens. All respiratory diagnoses were determined by the treating physician and were independent of study staff. Controls consisted of siblings of the cases or randomly selected children from the community. We defined “admission to hospital” as admission to a non-ED hospital ward. We also collected 19 convalescent/quiescent nasal samples from children who were followed-up at least 6 wk after an acute exacerbation of asthma or wheeze; only a subset of these samples (5/19) were paired with acute samples. The hospital’s human ethics committee approved the study (MAVRIC approval 1761 EP), and parental/guardian written informed consent was obtained prior to recruitment. An independent sample of children from within the MAVRIC cohort was used as a validation cohort.

A detailed questionnaire and medical records were used to provide demographic and medical information for each participant during recruitment and follow-up. A child was considered positive for aeroallergy if they had a positive response to either 1) the skin prick test completed on nine allergens (rye grass, mixed grasses, dog, cat, cockroach, Alternaria tenuis, Aspergillus fumigatus, Dermatophagoides farinae, Dermatophagoides pteronyssinus) or 2) a positive specific IgE to either cat or house dust mite (>0.35 kU/L) at either the acute or convalescent visit or positive answers to the acute questionnaire for the questions “Does your child suffer from hay fever?” or “Does your child suffer from allergies to: 1) grasses/pollens or 2) dust mite?” A child was considered positive for allergy overall if they had 1) aeroallergy, 2) a positive response to the skin prick tests for either cows’ milk or egg white, 3) a parental report of a history of anaphylaxis, or 4) a high total IgE at either the acute or convalescent visit.

Acute asthma severity scores were assigned to each child at recruitment using a modified National Institutes of Health score for children over 2 y of age (12) and included assessments of respiratory rate appropriate for age, oxygen saturation, auscultation, retractions, and dyspnea. A separate, age-appropriate severity score was used for children aged under 2 y (13) with assessments of heart rate, respiratory rate, wheezing, and accessory muscle use. Separate severity z-scores were calculated for each child within each of the two age groups to provide standardized scores across the whole cohort.

The time to next representation or admission to a public hospital in Western Australia with any acute respiratory illness for each child, within the first and second year of observation, was obtained from the state hospital database retrospectively.

A nasal secretion sample from each participant was collected to test for the presence of respiratory viruses and bacteria. Nasal swab specimens were obtained from each child using flocked swabs (Copan, Italy) and were taken immediately to the laboratory for processing for gene expression profiling.

Common respiratory pathogens (adenovirus; respiratory syncytial virus [RSV] types A and B; parainfluenza viruses 1–4; influenza viruses A, B, and C; metapneumovirus; Bordetella species; Mycoplasma pneumoniae; Chlamydophila pneumoniae; Haemophilus influenzae; Pneumocystis jirovecii; Staphylococcus aureus; Streptococcus pneumoniae; and pyogenes) were identified using a tandem multiplex real-time PCR assay as previously described (14). RV was detected and genotyped by a molecular method as previously described (15, 16).

Total RNA was extracted from nasal swabs using TRIzol (Invitrogen) followed by RNeasy MinElute (QIAGEN) (10). The quality and integrity of the RNA was analyzed on the nanodrop and bioanalyzer (Agilent). Total RNA samples were shipped on dry ice to the Ramaciotti Centre for Genomics at the University of New South Wales, Sydney, Australia, for processing and hybridization to Human Gene 2.1 ST microarrays (Affymetrix). The raw microarray data are available from the Gene Expression Omnibus repository (https://www.ncbi.nlm.nih.gov/geo; accession number: GSE103166).

The microarray data were preprocessed in R employing the Robust Multi-Array algorithm (17). A custom chip description file (hugene21sthsentrezg; version 20) was used to map probe sets to genes (18). The quality of the microarray data was assessed employing the R package arrayQualityMetrics. Low-quality/outlying arrays were identified on the basis of relative log expression and normalized unscaled standard error metrics and removed from the analysis (19, 20). Differentially expressed genes were identified using Linear Models for Microarray Data (LIMMA) in conjunction with surrogate variable analysis (SVA) (21, 22). LIMMA entails fitting a linear model to the log expression intensities for each gene, and Empirical Bayes methods are employed to moderate the SE estimates, resulting in improved power and more-stable inference (21). Surrogate variable analysis estimates and captures all sources of hidden biological and/or technical variation that may potentially confound the analysis, and the estimated surrogate variables are added as covariates to the LIMMA models. Where applicable, the correlation between paired samples was accounted for in the linear model fit through use of the duplicateCorrelation function. The p values derived from the LIMMA models were adjusted for multiple testing employing the false discovery rate method (21). Noninformative probe sets were identified using the proportion of variation accounted by the first principal component (PVAC) algorithm and filtered out of the results from the differential expression analysis (23).

Molecular phenotypes were identified using consensus hierarchical clustering (24). This algorithm employs data resampling techniques to find consensus across thousands of runs of a cluster analysis (24). Prior to cluster analysis, batch-like effects and other sources of unwanted variation were modeled and removed using the RUVnormalize algorithm (25). This algorithm leverages a set of negative control genes that are not associated with the outcome of interest to model unwanted variation and regress it out of the data. To select negative control genes, genes associated with wheezing exacerbations were first identified in case/control comparisons between wheezing/convalescent or wheezing/control samples using LIMMA/SVA, and those genes with an adjusted p value <0.2 for any comparison were excluded from control gene selection. Negative control genes (n = 1000) were selected from the remaining genes using the empNegativeControls function from the RUVCorr package (26). After removal of unwanted variation using the RUVnormalize algorithm, we demonstrated that principal components of gene expression were no longer strongly correlated with estimated proportions of specific cell types. The RUVnormalized gene expression data were converted to z-scores, and consensus hierarchical clustering was performed using the set of genes that was associated with the outcome of interest (adjusted p value <0.1). Hierarchical cluster analysis was performed using Pearson correlation and ward linkage.

Pathways analysis was performed using Enrichr software and the Reactome database (27).

Gene networks were constructed employing experimentally supported findings from published studies curated in the Ingenuity Systems Knowledge Base (Ingenuity Systems, Redwood City, CA). Direct and indirect molecular relationships were considered from all available categories, including activation, inhibition, localization, modification, molecular cleavage, phosphorylation, protein–DNA interaction, protein–protein interaction, regulation of binding, transcription, translation, and ubiquitination. The circular layout was employed to display the network graph object, and hubs were positioned at the center of the network. The nodes were colored based on the expression ratio; red nodes denote upregulated genes, and green nodes denote downregulation.

Upstream regulator analysis (Ingenuity Systems) was employed to infer the molecular drivers of the observed differential gene expression patterns (28). Two statistical metrics were calculated. The overlap p value measures enrichment of known downstream target genes for a given upstream regulator among the differentially expressed genes. The activation z-score is a measure of the pattern match between the direction of the observed gene expression changes (in terms of up/downregulation) and the predicted pattern based on prior experimental evidence. Pathways with an activation z-score >2.0 are predicted to be activated, and an activation z-score < −2.0 indicates pathway inhibition. If the activation z-score lies between −2.0 and 2.0, the activation state of the pathway cannot be predicted.

Cell type enrichment was examined in microarray-based gene expression profiles using the xCell webtool (University of California, San Francisco) (29). Gene signatures for 64 cell types were determined from >1800 transcriptomic profiles of purified cell samples, allowing reliable enrichment analysis to investigate tissue heterogeneity. Cell types were chosen for further analysis based on their potential importance in the upper airways. Enrichment scores were converted to proportional cellular composition, allowed for by the linearity assumption, so that the sum of the proportions was equal to one. Statistical comparisons were performed using nonparametric Kruskal-Wallis and Mann–Whitney U tests (GraphPad software; La Jolla, CA).

Total RNA extracted from nasal swab samples was reverse transcribed into cDNA using the QuantiTect Reverse Transcription Kit (QIAGEN). qPCR was performed using QuantiFast SYBR Green PCR Master Mix (QIAGEN) on the ABI 7900HT Sequence Detection System (Life Technologies). Primer sequences were obtained from Primerbank (30) and purchased from Sigma-Aldrich. Standard curves were prepared from serially diluted quantitative RT-PCR (qRT-PCR) products. Melt-curve analysis was conducted for all samples to confirm the specificity of amplified products. Relative expression levels of target genes were calculated by normalization to the housekeeping genes HMBS, PPIA, and PPIB (9). The qRT-PCR was log2 transformed, converted to z-scores, and analyzed by consensus hierarchical clustering.

Comparisons of all categorical variables between the clusters were performed using χ² tests, whereas ANOVA was used for continuous variables. Kaplan–Meier survival curve was used to assess the time to first representation or admission after recruitment. A p value <0.05 was considered significant. All analyses were performed using SPSS version 22 (Chicago, IL).

The study was based on a case/control design. The cases (n = 56) consisted of children who presented to the ED with an acute exacerbation of asthma or wheeze. The controls consisted of children who were either siblings of the cases or they were recruited from the general community (controls, n = 31). Convalescent samples (n = 19) were available from children after they had recovered from an acute exacerbation of asthma or wheeze, but only a subset of these convalescent samples (n = 5/19) were paired with acute samples. Samples from an independent group of children (n = 99) with exacerbations of asthma or wheeze were used as a validation cohort. The characteristics of the study participants are presented in Table I.

Children with wheezing exacerbations were younger (mean = 4.35 [SD: 3.31] y) than controls (mean = 6.40 [SD 4.41] y, p = 0.025) and fewer were white (46.4 versus 78.6%, p = 0.006). A higher proportion of the children with wheezing exacerbations were sampled during winter, compared with convalescence (p = 0.001) and controls (p = 0.028). Respiratory virus, in particular, RV, was widely detected during wheezing exacerbations (83.9 and 66.1%, respectively) compared with convalescence (31.6 and 26.3%, p < 0.001 and p = 0.003, respectively) and controls (45.2 and 32.3%, p < 0.001 and p = 0.003, respectively). RV species C was more prevalent during wheezing exacerbations (46.4%) compared with controls (12.9%, p = 0.002). The number of viruses detected during acute exacerbations (mean = 0.96 [SD 0.60]) was higher than during convalescence (mean = 0.37 [SD 0.60], p < 0.001). No difference in bacterial detection was observed between the groups.

The proportion of children with wheezing exacerbations recruited in winter in the validation cohort was lower than the discovery cohort (p < 0.001, Table I). The validation cohort was also more atopic (78.8 versus 60.7%, p = 0.024) and had lower detection rates of respiratory viral and bacterial pathogens compared with the latter. Respiratory symptoms, hospital admissions, and medication usage were not different in the discovery and validation cohorts.

Nasal swab specimens were collected from the children, and gene expression patterns were profiled on microarrays. Cellular composition data of the samples were not available, and therefore a computational approach was employed to estimate the proportions of different cell types directly from the microarray profiles (31). As illustrated in Fig. 1A, the highest proportions were observed for neutrophils, epithelial cells, and monocytes. Stratification of the subjects by case/control status and RV detection revealed that wheezing exacerbations were associated with increased proportions of M1 macrophages and decreased proportions of conventional dendritic cells (Fig. 1B). The proportions of the other cell types were mostly comparable across the groups. The relationship between case/control status, RV detection, and cellular composition for each individual child is illustrated in Fig. 1C.

To identify differentially expressed genes, we undertook group-wise comparisons employing LIMMA in combination with surrogate variable analysis, which adjusts the analysis for all estimated sources of unwanted biological and technical variation (see 2Materials and Methods). Comparison of gene expression patterns between children with RV-positive wheezing exacerbations (n = 37) versus RV-negative controls (n = 21) revealed that 137 genes were upregulated and 82 genes were downregulated (adjusted p value <0.05, Supplemental Table I). Pathways analysis demonstrated that the upregulated genes were highly enriched for type I IFN signaling (adjusted p value = 4.8 × 10⁻²², relevant genes are highlighted in Supplemental Table I). Type I IFN signaling (adjusted p value = 8.5 × 10⁻³¹, Supplemental Table II) was also upregulated in children with RV-positive exacerbations in comparison with RV-negative convalescent children (n = 14).

Gene expression patterns were then assessed in children with RV-negative exacerbations. Comparison of these children (n = 19) with RV-negative controls (n = 21) revealed that 17 genes were upregulated and 10 genes were downregulated (adjusted p value <0.05, Supplemental Table III). Similar results were obtained by comparing children with RV-negative exacerbations versus RV-negative convalescent children (n = 14), in which 30 differentially expressed genes were identified (adjusted p value <0.05, Supplemental Table IV). Although a subset of these genes was also upregulated in children with RV-positive exacerbations (e.g., IL-18R1, CD163), a notable difference was the lack of an IFN signature.

The findings above suggested that children with RV-positive exacerbations mount an IFN response, but there appears to be additional heterogeneity. As the data analysis strategy relied on group-wise comparisons (with LIMMA/SVA), this approach cannot shed any insight into patient-to-patient variations in gene expression patterns within each group. To obtain more detailed information in this regard, we employed hierarchical clustering (see 2Materials and Methods). To ensure that the clustering did not simply reflect variations in cellular composition (observed in Fig. 1), we used a set of negative control genes (i.e., genes not related to the outcome of interest) to model unwanted variation in the data and removed these effects using regression (see 2Materials and Methods and Supplemental Fig. 1). This strategy preserves the gene expression signature associated with acute exacerbations in the data and removes the strong correlation structure between samples that reflect variations in cellular composition (Supplemental Fig. 1). As illustrated in Fig. 2A, cluster analysis of the corrected data divided the subjects into five distinct clusters or molecular phenotypes (labeled S1–S5). Likewise, the genes were also divided into five clusters labeled G1–G5. Notably, the majority of children with exacerbations (45/56 = 80%) were found within clusters S1 and S2, and there was no difference in RV detection between these cluster groups of children (RV detection rates for children with exacerbations in cluster S1 versus cluster S2 were 73 versus 58%, p value = 0.347). The other three clusters, S3–S5, contained the bulk of the controls and the convalescent samples as well as 11 remaining children with exacerbations.

To examine the overall expression of each gene cluster across the children, principal component analysis was employed to summarize the expression data for each cluster, and the first principal component was plotted across the subjects, stratified into their respective clusters. This analysis revealed that the most striking difference between the subjects in cluster groups S1 and S2 was the differential expression of the set of genes in cluster G5 (Fig. 2B), and this cluster of genes was strongly enriched for type I IFN signaling (data not shown). Moreover, prior knowledge–based reconstruction of the wiring diagram of the underlying gene networks for each cluster revealed that IRF7, a master regulator of type I IFN responses (32), was the dominant hub gene identified in cluster G5. Accordingly, we designated the children in clusters S1 and S2 as IRF7^hi versus IRF7^lo molecular phenotypes.

To further characterize the IRF7 phenotypes, children with exacerbations were stratified into IRF7^hi (n = 26) and IRF7^lo (n = 19) subgroups and compared with RV-negative controls employing LIMMA/SVA. We followed this strategy because a direct comparison of IRF7^hi/low phenotypes would only reveal differences between the respective responses, and we wanted to know if the phenotypes operate through discrete and/or overlapping pathways. The data showed that 208 genes were upregulated and 157 genes were downregulated in children with IRF7^hi exacerbations compared with controls (Fig. 3A, left panel; Supplemental Table V). IRF7^lo exacerbations were characterized by upregulation of 96 genes and downregulation of 31 genes (Fig. 3A, right panel; Supplemental Table VI). These analyses revealed an overlapping response signature, which comprised 52 upregulated genes and 11 downregulated genes. Gene network analysis revealed that IRF7 gene networks were upregulated in the children with IRF7^hi exacerbations (Fig. 3B, left panel). IRF7^lo exacerbations lacked an IRF7 signature and instead were characterized by upregulation of Th2-associated pathways (e.g., IL-4R, FCER1G, ARG1) and downregulation of IFN-γ (Fig. 3B, right panel). We also built a combined network to illustrate the unique and overlapping gene network patterns for each phenotype (Fig. 4).

Next we employed upstream regulator analysis to identify molecular drivers of IRF7^hi and IRF7^lo exacerbation responses (28, 33). The data showed that IRF7^hi exacerbations were putatively driven by upregulation of IFNL1, IRF7, and IFN-α signaling (Fig. 3C, left panel, Supplemental Table VII). In contrast, IRF7^lo exacerbations were predicted to be driven by upregulation of cytokine and growth factor signaling pathways (e.g., IL-6, IL-10, TGF-β, ERK, CSF3, EGF; Fig. 3C, right panel, Supplemental Table VIII). It is noteworthy that multiple inflammatory pathways (e.g., IL-1β, IL-2, IL-4, IL-13, TNF, OSM, IFN-γ) were featured in the upstream regulator analysis results for IRF7^lo exacerbations; however, the activation z-scores were not significant, indicating that the direction of the gene expression changes in terms of up- and downregulation were not consistent with the known role of these pathways in the regulation of gene expression (28).

We reexamined the cellular composition data from Fig. 1 in IRF7^hi and IRF7^lo exacerbation responses, and we found that both phenotypes were associated with decreased proportions of conventional dendritic cells and similar proportions of M1 macrophages (Supplemental Fig. 2). Stratification of IRF7 phenotypes by RV detection revealed further heterogeneity, including upregulated proportions of Th2 cells in children with RV-positive, IRF7^lo exacerbations (Supplemental Fig. 3).

To identify candidate pathways that potentially link exacerbation responses with expression of asthma-related traits, we employed Pubmatrix (34) to screen the literature for relevant studies (Supplemental Table IX). We found that multiple genes that were upregulated in children with IRF7^hi exacerbations (e.g., CASP1, CD274, CXCL11, DDX58, haptoglobin, IFIH1, IRF7, P2RX7, PHF11, selectin L, SERPING1, STAT1, TLR4, TLR5, TNFAIP3, TNFAIP6, TNFSF13B) or IRF7^lo exacerbations (e.g., AGR3, C3AR1, EPCAM, IL4R, LDLR, NCR1, OCLN, SERPINB2, SERPINB3, THBS1, TIMP1, VCAN) have been previously studied in the context of asthma and/or related traits. Moreover, overlap signature of genes that were commonly upregulated in children with IRF7^hi and IRF7^lo exacerbations has also been previously studied in this context (e.g., ADAM17, ARG1, ARG2, CD163, CD38, FCER1G, IL18R1, PLA2G4A, S100A12, TLR2).

We next investigated the biological and clinical characteristics of IRF7^hi and IRF7^lo exacerbations. These groups were not related to any technical variables measured in the study (Supplemental Table X). There was also no difference in season of recruitment, the detection of viral or bacterial pathogens, or the use of medications (Table II). The prevalence of atopy, including allergy to aeroallergens, was higher in IRF7^hi exacerbations (76.9 and 73.1%, respectively) compared with IRF7^lo exacerbations (42.1 and 31.6%, p = 0.029 and 0.008, respectively). Children with IRF7^lo exacerbations presented to the hospital much later after initial symptoms (mean = 4.74 [SD 4.03] d) than IRF7^hi exacerbations (mean = 2.31 [SD 1.98] d, p = 0.011). This was also reflected in the duration of cough prior to hospitalization (IRF7^lo: 5.62 [SD 3.20] d and IRF7^hi: 1.96 [SD 1.51] d, p = 0.000027). Children with IRF7^lo exacerbations were at least 4.6 times more likely to be admitted to the hospital compared with IRF7^hi exacerbations (odds ratio: 4.65, p = 0.018). Time to subsequent first representation/admission to the hospital with an exacerbation was shorter in IRF7^lo exacerbations compared with IRF7^hi exacerbations (Supplemental Fig. 4, Table II). Within the first year of follow-up, more children with IRF7^lo (68.4%) represented/readmitted to the hospital with a respiratory exacerbation compared with those with IRF7^hi (30.8%, p = 0.017). All associations remained significant after adjustment for age, gender, aeroallergen allergy, and white ethnicity.

To determine if IRF7^hi and IRF7^lo exacerbation phenotypes and their clinical correlates could be found in another group of children, qRT-PCR was employed to measure gene expression patterns in nasal swab samples from an independent sample of children. To select a panel of genes for qRT-PCR analysis, we focused on genes that were representative of IRF7^hi exacerbations (DDX60, IFNL1, IRF7, ISG15, Mx1, RSAD2, and the downregulated gene IL-33; Fig. 4), IRF7^lo exacerbations (NCR1, THBS1; Fig. 4), or that were common to both phenotypes (ARG1, CD163, FCER1G, IL-1R2, IL-18R1, TLR2; Fig. 4). We also demonstrated that cluster analysis of this restricted panel of genes could largely reproduce the clustering result based on the full microarray data set (Supplemental Fig. 5). The qRT-PCR data were normalized to three endogenous control genes (HMBS, PPIA, PPIB), creating three separate variables for each gene, which were used for consensus clustering. The analysis segregated the subjects into five clusters (Fig. 5); three of these had elevated expression of IRF7/IFN responses (combined as IRF7^hi; black dendrogram in Fig. 5), and the remaining two clusters had low IRF7/IFN responses (IRF7^low1 and IRF7^low2; red and green dendrograms, respectively, in Fig. 5). A longer time lag was observed from first symptoms to hospital presentation in the IRF7^low1 subjects (mean = 3.90 [SD 4.59] d) compared with the IRF7^hi group (mean = 2.20 [SD 1.76] d, p = 0.023). These symptoms included cough (p = 0.015), wheeze (p = 0.022), and shortness of breath (p = 0.02). Fever was more prevalent in the IRF7^hi group (59.2%) compared with the IRF7^low1 and IRF7^low2 groups (24.1%, p = 0.004 and 10.0%, p < 0.001, respectively). Runny nose was also more prevalent in the IRF7^hi group (75.5%) compared with the IRF7^low2 group (45.0%, p = 0.024). All associations remained significant after adjusting for age, gender, aeroallergy, and white ethnicity.

Our study is the first, to our knowledge, to investigate acute wheeze/asthma exacerbation phenotypes in children using a systems biology approach. Our computational analyses identify IRF7 as a gene network hub and reveal two distinct IRF7 molecular phenotypes: IRF7^hi exhibiting robust upregulation of the Th1/type I IFN response and IRF7^lo with an alternative activation signature marked by upregulation of cytokine and growth factor signaling and downregulation of IFN-γ. Importantly, the two phenotypes also produced distinct clinical phenotypes. Compared with children with IRF7^hi, those with IRF7^lo had a slower progression of illness from initial symptoms to hospital presentation, a greater likelihood of a hospital admission, and a greater chance of representation with a further exacerbation. Exacerbation severity at presentation was not different between the two IRF7 patterns, perhaps because presentation to the hospital is likely to be determined by symptoms reaching a severity threshold that causes parental concern. However, also worth noting is that the apparently impaired IRF7 response in IRF7^lo children was not associated with an increase in exacerbation severity at presentation, although perhaps expectedly the reduced IRF7 response was associated with slower resolution of the episode. These findings were unaffected by the respiratory viruses or bacteria detected and medication use. Investigation of a validation cohort using qRT-PCR produced similar findings despite variations in the patient characteristics between the discovery and validation cohorts. In summary, our findings reveal two distinct phenotypes with clear differences in gene regulation patterns, either upregulation of robust innate immune responses or cytokine and growth factor signaling, and associated differences in clinical characteristics.

Asthma and wheezing exacerbations are largely triggered by RV infections, but the underlying mechanisms are not well understood (2). Previous in vitro studies suggested that RV-induced IFN responses were deficient in bronchial epithelial cells from subjects with asthma, resulting in increased viral loads and exaggerated secondary responses (3, 4, 35). However, in vivo studies of immune response patterns in the airways of both children and adults with RV-induced exacerbations found that IFN responses were increased, not deficient (9–11, 36). Moreover, a clinical trial that evaluated the use of inhaled IFN-β therapy in adults at the first signs of a cold to prevent worsening of their asthma symptoms failed to achieve its primary end point (37). Our data extend these previous findings by demonstrating that wheezing exacerbations in children are heterogeneous and characterized by two very different IRF7 molecular phenotypes. Upstream regulator analysis suggested that IRF7^hi exacerbations were driven by upregulation of IFNL1, IRF7, and IFN-α signaling. In contrast, IRF7^lo exacerbations were putatively driven by upregulation of cytokine and growth factor signaling (i.e., IL-6, IL-10, TGF-β, CSF3, EGF) and downregulation of IFN-γ.

IRF7 is a master regulator of type I and type III IFN gene expression (32, 38–40). We previously reported that IRF7 gene networks were upregulated in nasal wash samples from asthmatic children experiencing mild-to-moderate viral exacerbations (10). We have also shown that IRF7 promotes RV-induced innate antiviral responses and limits IL-33 mRNA expression in cultured bronchial epithelial cells (33), and furthermore, in our current study, IL-33 was downregulated in IRF7^hi exacerbation responses. Girkin et al. (40) reported that knockdown of IRF7 abolished RV-induced type I IFN responses in the airways in a mouse model. Werder et al. (41) reported that anti–IL-33 therapy suppressed type 2 inflammation, boosted antiviral immunity, and decreased viral replication in a mouse model of virus-induced asthma exacerbations. Together, these findings highlight a mutually antagonistic role for IRF7-related antiviral immunity and IL-33–driven type 2 inflammation in driving alternative exacerbation phenotypes.

Children with IRF7^lo exacerbations in the discovery cohort exhibited a delayed progression from first symptom onset to hospital presentation. Given that our study design entailed recruitment of children at ED presentation, we could not determine if IRF7 gene networks were initially upregulated closer to the onset of infection and subsequently waned or, alternatively, if they were never upregulated in the first place. To address this issue, an alternative study design would be required that entails regular sampling of exacerbation-prone children during the RV season (42). It would also be of interest to further study the stability of these identified phenotypes to learn if subjects experience the same type of response over multiple wheezing exacerbation events or not. Notwithstanding this, our validation cohort comprised approximately twice as many cases as the discovery cohort, and this larger sample enabled the identification and characterization of two subgroups of IRF7^lo children, and only one of these subgroups was characterized by delayed progression. The activation of growth factor signaling pathways and downregulation of IFN-γ may reflect the immune response entering a resolution phase, however, at the same time these children were symptomatic. Moreover, many were also RV positive, and given that TGF-β signaling promotes RV replication, upregulation of this pathway may prolong infection and delay viral clearance (43, 44). It is also known that frequent severe exacerbations are associated with deficits in lung function growth (children) and accelerated lung function decline (adults) (45). Thus, repeated cycles of inflammation, growth factor signaling, and repair may alter the structure and function of the airways underlying this phenotype.

The mechanisms that determine expression of disease symptoms among children with IRF7^hi and IRF7^lo exacerbations are unknown. Given that children with IRF7^hi exacerbations elicit robust antiviral responses, one possibility is that the airways of these children are sensitive to the host response to respiratory viruses. In this regard, although we acknowledge that the expression profiling studies were performed in the upper airways and asthma symptoms manifest in the lower airways, several pathways associated with the IRF7^hi phenotype are known to impact respiratory function. PD-L1 (encoded by CD274) is an immune checkpoint that delivers an inhibitory signal for T cell activation. Upregulation of PD-L1 during respiratory bacterial infections in early life suppresses the IL-13 decoy receptor IL-13Ra2, resulting in persistent airway hyperresponsiveness (46). MDA5 (encoded by IFIH1) is a pattern recognition receptor that senses RV-derived RNA. MDA5-deficient mice infected with RV have delayed type I IFN responses, impaired type III IFN responses, and reduced airway hyperresponsiveness (47). The proinflammatory effectors TNF, IFN-γ, and IFN-γ–induced protein 10 can also promote airway hyperresponsiveness in animal models (48–50).

Another possibility is that IRF7^hi and IRF7^lo exacerbation responses converge on a final common pathway to precipitate respiratory symptoms (Fig. 4). For example, CD38 is a receptor with enzymatic activity, which hydrolyses NAD, generating reaction products that modulate calcium signaling. It is expressed on immune and airway smooth muscle cells, and it plays a dual role in asthma by enhancing airway inflammation and contractile responses in smooth muscle (51). FCER1G is a component of the high-affinity IgE receptor. Anti-IgE therapy neutralizes serum IgE, reduces the expression of the high-affinity IgE receptor on dendritic cells and mast cells, and it also reduces the frequency of asthma exacerbations (52, 53). Phospholipases A2 are involved in the generation of eicosanoids from arachidonic acid. Knock-in of human sPLA2 into mice enhances airway inflammation and airway hyperresponsiveness (54). TLR2 is a pattern recognition receptor that acts as a sensor for RV capsid (55). In a mouse model of combined RV and allergen exposure, TLR2 signaling in macrophages was required for the induction of airway inflammation and airway hyperresponsiveness (56). In summary, our data have identified multiple candidate pathways that potentially link IRF7^hi and IRF7^lo exacerbation responses with expression of respiratory symptoms, and these pathways represent logical candidates for future drug development programs.

This study has limitations that should be acknowledged. The expression profiles were derived from nasal swab samples, and it is not known to what extent the data reflect inflammatory mechanisms operating in the lower airways. The nasal swab samples comprised a mixed cell population. Follow-up studies employing focused analyses in isolated cell types or single-cell transcriptomics will be required to dissect the role of specific subpopulations of cells in this disease. The study participants were sampled during natural exacerbations, and it is not possible to control for all of the variables that may potentially impact the data (e.g., age, gender, ethnicity, natural allergen exposure, pathogen strains and combinations, medications). To address this issue, our analysis strategy employed surrogate variable analysis to systematically estimate and adjust the analysis for all sources of hidden biological and/or technical variation. Some of the clinical characteristics associated with IRF7 phenotypes in the discovery cohort did not replicate in the validation cohort. This may be due in part to variations in the demographics and clinical characteristics between the two cohorts as well as the fact that IRF7 phenotypes were defined in the discovery cohort by microarray analysis of a large number of genes, whereas in the validation cohort it was based on qRT-PCR analysis of a restricted gene panel. The use of microarray technology for expression profiling is also a limitation because RNA sequencing can detect a lot more transcripts (57). Finally, although our analyses have characterized gene network patterns underlying exacerbation phenotypes and unveiled candidate molecular drivers of the responses, further studies will be required to dissect the mechanisms that give rise to these phenotypes and drive the expression of respiratory symptoms. Notwithstanding these limitations, our findings demonstrate that exacerbation responses in children are heterogeneous and comprise IRF7^hi versus IRF7^lo molecular phenotypes that determine clinical phenotypes. Future clinical trials targeting the IFN system in this disease should be stratified on the basis of these IRF7 phenotypes.

We thank all the children and families who agreed to take part in the study.

The authors have no financial conflicts of interest.