Decreased activation of inflammatory networks during acute asthma exacerbations is associated with chronic airflow obstruction.

Asthma exacerbations are associated with subsequent deficits in lung function. Here, we tested the hypothesis that a specific pattern of inflammatory responses during acute exacerbations may be associated with chronic airway obstruction. Gene coexpression networks were characterized in induced sputum obtained during an acute exacerbation, from asthmatic children with or without chronic airflow limitation. The data showed that activation of Th1-like/cytotoxic and interferon signaling pathways during acute exacerbations was decreased in asthmatic children with deficits in baseline lung function. These associations were independent of the identification of picornaviruses in nasal secretions or the use of medications at the time of the exacerbation. Th2-related pathways were also detected in the responses, but variations in these pathways were not related to chronic airways obstruction. Our findings show that decreased activation of Th1-like/cytotoxic and interferon pathways is a hallmark of acute exacerbation responses in asthmatic children with evidence of chronic airways obstruction.

Introduction

In most cases of asthma, the first symptoms of the disease occur during childhood,1 and in a subgroup of children, the disease is associated with deficits in airway function that track with age and predict persistence of asthma symptoms into adult life.2–4 There is now strong evidence suggesting that these deficits are acquired during the course of the disease,5 are more likely to become apparent during childhood,6 and that chronic administration of inhaled corticosteroids does not prevent them.7 Moreover, at least one-third of older asthmatics develop chronic obstructive pulmonary disease (COPD),8 and among that third, airflow limitation is often present in early adult life.9 Thus, a better understanding of the molecular mechanisms that may be associated with airflow limitation in children with asthma could help design new treatments to prevent persistent asthma and asthma-associated COPD.

Recent longitudinal studies have suggested that, among children and adults with asthma, those with a higher incidence of acute exacerbations are more likely to have subsequent deficits of lung function growth10 or excess decline,11 respectively. The mechanisms that underlie the associations between exacerbations and chronic lung disease in humans are unknown. Evidence of infection with rhinoviruses has been reported in 50–70 % of asthma exacerbations,12 and accumulating data suggests that immune responses to viruses may be altered in asthmatics. When infected with rhinovirus, bronchial epithelial cells obtained from adult atopic asthmatics have deficits in interferon type I and III responses.13, 14 Impaired interferon gamma and augmented Th2 responses in blood and BAL T cells from asthmatics is associated with rhinovirus-induced clinical illness severity and viral load.15 Rhinovirus-induced interferon gamma-to-IL-5 response ratios in PBMC cultures in vitro are positively correlated with forced expiratory volume in 1 second (FEV1),16 and these same response ratios in sputum are inversely related to viral clearance in vivo.17 A recent study in mice demonstrated that respiratory viral infections can trigger the development of a chronic asthma/COPD-like lung disease after the virus is cleared to trace levels.18 In this mouse model, invariant natural killer T (iNKT) cells play a central role in disease pathogenesis; nevertheless, there are conflicting data on the role of iNKT cells in the pathogenesis of human asthma.19, 20

We postulated that children whose asthma is associated with chronic airflow limitation could have a specific pattern of inflammatory responses in the airways that predisposes them to the development of deficits in airway function.21 To test this hypothesis, we induced sputum at the time of a moderate exacerbation and 7–14 days later in a group of children with asthma followed prospectively, and compared global patterns of gene expression in sputum cells in those with or without baseline airflow limitation. Our results show that the activation of Th1-like/cytotoxic and interferon signalling pathways during exacerbations is decreased in asthmatic children with evidence of chronic airflow obstruction. Moreover, regardless of the presence or absence of chronic airflow obstruction, acute exacerbations were associated with markers of iNKT cells, and the latter were highly correlated with cytokines that promote Th1/cytotoxic responses (IL-12A, IL-21).

Methods

Study population and protocol design

The study population consisted of 218 mild/moderate persistent asthmatic children aged 6–18 years who had at least one acute asthma exacerbation during the previous year. The study protocol included three visits. At enrollment, subjects had to be symptom-free; asthma severity and control was assessed by a study physician and adjustments were made, if necessary, to achieve national guideline recommendations for control. Baseline spirometry was assessed using an appropriately calibrated Jaeger spirometer (Jaeger, Erich, Jaeger, Wurzburg, Germany). The ratio between FEV1 and forced vital capacity (FVC) was used to assess airflow limitation at baseline, as suggested by Rasmussen et al.22 Skin test reactivity to local allergens was measured (Alternaria, dust mite mix [Dermatophagoides farinae plus Dermatophagoides pteronyssinus], olive, Bermuda, careless weed, mulberry, cockroach, cat, dog, mouse, and ragweed). After enrollment, the children were followed for 18 months or until their first asthma exacerbation, whichever came first.

A moderate exacerbation was defined with criteria equivalent to those recently proposed by the ATS/ERS consensus.23 If a child experienced symptoms of cough, dyspnea, chest tightness, and/or wheeze, they were instructed to initiate use of albuterol (2 puffs, 90 mcg/puff) by MDI every 20 minutes for up to 1 hour and then every 4 hours if necessary. For those who could use peak flow meter, a reading of < 80 % of personal best was considered indicative of an exacerbation. If the subject could not get symptom relief after 3 treatments, parents were instructed to contact the study center immediately. These three conditions (symptoms of exacerbation, low peak flow readings, or lack of relief or persistence of symptoms after three treatments) were considered indicative of an exacerbation. Participants who met these criteria were scheduled for a visit at the study center within 24 hours. Ascertainment of an exacerbation was ultimately made by a study physician. This research was approved by the Institutional Review Board of the University of Arizona and informed consent was obtained for all participants.

Sputum induction

At the time of the acute exacerbation visit, a physical examination was performed by a study physician, lung function was assessed by spirometry, nasal secretions were collected to test for evidence of a picornavirus infection (online methods), and sputum was induced based on the techniques recommended by Gershman et al24 with slight modification. Following the measurement of spirometry pre- and post-bronchodilator, a peak flow measurement was made to determine the subject’s baseline. Participants then inhaled 3 % saline for 1.5 minutes, after which any accumulated saliva was discarded, followed immediately by coughing and expectoration of sputum. Peak flow was checked to assure that there had not been a 15 % or greater drop following saline inhalation. This procedure was repeated at 2 minute intervals six times. Following the final inhalation of 3 % saline and expectoration, a full spirometric maneuver was done in place of peak flow to assure no decline from baseline values. Only children whose FEV1 was > 70 % of predicted were eligible for sputum induction. Seven to fourteen days after the acute episode, children were seen again in the study clinic and physical examination, nasal washes, and sputum collection were repeated as described above.

Sputum processing and RNA stabilization

Sputum was stored at 4 °C and processed within one hour of collection. Briefly, whole sputum was diluted with a volume of 0.1 % Dithiothreitol (DTT-Sputolysin 10 %; Calbiochem Corp, La Jolla, CA) equal to that of the weight of the original sample in grams, and incubated on a shaking water bath at 37 °C for 15 minutes with intermittent aspiration via pipetting every 5 minutes. Homogenized sputum was centrifuged (800 g) for 10 minutes. The cell pellet was immediately resuspended in RNAlater Stabilization Reagent (Qiagen, Valencia, CA) as per manufacturer’s recommendations. Total cell counts and differentials were determined from an aliquot of the original homogenized sample. Slides were prepared (Cytospin; Shandon; Runcorn, UK) and stained with a Wright-Giemsa stain. Differential cell counts were made by a blinded observer. One hundred cells were counted for each sample. Differential cell counts are expressed as percentages of total cells.

Microarray-based expression profiling studies

Total RNA from sputum samples stored in RNALater was extracted with TRIzol (Invitrogen, Carlsbad, CA) followed by RNeasy (QIAgen, Valencia, CA). In preliminary studies based on Bioanalyzer analysis, we noted some variation in the integrity of RNA from sputum. However, the microarray and real time quantitative reverse transcription PCR (qRT-PCR) protocols employed in the study are based on random priming and are thus tolerant to variations in RNA quality. Total RNA samples (n=20) were labelled and hybridized to Human Gene ST 1.0 microarrays (Affymetrix, Santa Clara, CA), at the Arizona Cancer Center Genomics Core, the University of Arizona. The microarray data are available from the Gene Expression Omnibus repository (accession GSE19903).

The microarray data was preprocessed in Expression Consol software (Affymetrix, Santa Clara, CA) employing the probe logarithmic intensity error algorithm with gc background subtraction, quantile normalization and iterPLIER summarization. The preprocessed data was imported in the R language for statistical computing (http://www.r-project.org/), and variance stabilization was performed by adding the small constant 16 to all the data points, followed by log2 transformation.

Reverse engineering gene network analysis

The microarray data was filtered to select highly variable genes (top 10 % on microarray – 3247 genes) as well as the top 1500 genes that differed in the respective responses (ie. in subjects with (n=10) or without (n=10) deficits in enrollment FEV1/FVC ratios) according to their statistical ranking from a Bayesian t.test analysis.25 Network analysis was then performed on the filtered dataset in all subjects employing the weighted gene coexpression network analysis (WGCNA) algorithm.26–28 The algorithm calculates absolute Pearson correlations for all pairwise gene-gene combinations across the test samples. The correlations are then raised to a power to emphasize stronger over weaker correlations. Genes that had a low overall correlation with the coexpression network were removed from the analysis (approx 25 % of initial genes removed). The topological overlap of the gene-gene correlations was calculated to quantify the extent that genes have similar overall patterns of correlations with other genes. The topological overlap similarity measure was subtracted from one to convert it into a distance measure and then analyzed by hierarchical clustering to group highly correlated genes into subnets (modules). The modules were defined from the dendrogram output of the cluster analysis employing an automated algorithm (cutreeDynamic). The overall expression of the modules was compared in the respective responses employing Gene Set Analysis without correcting for multiple testing.29 To determine if the modules were stable, a randomly selected sample was removed from the analysis just prior to the hierarchical clustering stage, and new modules were defined. This process was repeated an additional four times, and the stability of the modules was calculated as the proportion of genes from the original cluster that were detected in the same cluster, averaged over the five iterations.

Bioinformatics analysis of molecular signatures

The list of genes in the Th1-like/cytotoxic pathway (Table S2) was interrogated for significant overlaps with the collection of 1,892 curated molecular signatures from the Molecular Signatures Database. The database contains annotated pathways from online databases, and molecular signatures from published microarray studies, thus captures a broad range of biological, cellular, and clinical states. Statistically significant overlaps were identified based on the Hypergeometric distribution.30 This analysis was performed online (http://www.broadinstitute.org/gsea/index.jsp).

qRT-PCR validation studies

Total RNA was reverse transcribed with a combination of random nonamers and oligo-dT priming using the Quantitect reverse transcription kit with integrated genomic DNA removal (QIAgen, Valencia, CA). qRT-PCR analysis was performed with Quantitect SyBr green (QIAgen, Valencia, CA) on the 7900 thermocycler (Applied Biosystems, Foster City, CA). The primer assay sequences for FCER1A were obtained from Primer Bank (http://pga.mgh.harvard.edu/primerbank/index.html), and all other assays were obtained from QIAgen. Quantification was based on the relative standard curve method, and standards were prepared for each assay by serial diluting qRT-PCR products. The specificity of the qRT-PCR assays was confirmed by dissociation curve analysis and by testing negative RT control reactions. To select a housekeeping gene for normalization of the qRT-PCR data, we interrogated the microarray data for the following set of house keeping genes: ACTB, ALAS1, B2M, EEF1A1, GAPDH, GUSB, HMBS, HPRT1, PGK1, PPIA, RPL13A, RPL27A, RPL37A, RPLP0, RRN18S, SDHA, TBP, TFRC, TUBB, YWHAZ. We selected HMBS31 because it had the lowest variance. We also confirmed by qRT-PCR that the variance of HMBS expression was an order of magnitude lower than that of one of the most stably expressed genes in the genome (i.e., EEF1A1)32 and several orders of magnitude lower than RRN18S (18S rRNA). It is noteworthy that HMBS is expressed at low levels,31 thus providing more adequate adjustments for variations in RNA quality than highly expressed house keeping genes.

Real-time qRT-PCR detection of iNKT cells and T cells in sputum

T cells and iNKT cells can be detected in sputum via qRT-PCR analysis of T cell receptor transcripts.20 Most human iNKT cells express a Vα24-Jα18 (encoded by TRAV10 and TRAJ18; note Jα18 was formally JαQ) T cell receptor alpha chain rearrangement paired with a Vβ11 beta chain (encoded by TRBV25-1).33 T cells including iNKT cells express the T cell receptor constant chain (TRBC2).20 The primer assay for detection of Vα24 was Vα24-F-AAGCATCTGACGACCTTCTTG, and Vα24-R-AACAGGACCTCTCCCAGTATC.20 Primer assays for Vβ11 (F – CCTCTGCTACGTGGGCTTT; R – GCCCATGGTTTGAGAACATT) and TRBC2 (F – AACCACTTCCGCTGTCAAGT; R – GCTGGTAAGACTCGGAGGTG) were designed with Primer3Plus (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi) to amplify ensemble transcript sequences (ENST00000390398, ENST00000390420), and concordance of ensemble sequences with published sequences from iNKT cells33 was confirmed using BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi).

Statistical design and methodology

This study was designed with a two-stage analytical strategy. First, a case-control approach was employed; two groups of 10 subjects each from both extreme tails of the FEV1/FVC distribution at enrollment were selected for the microarray profiling studies. In the second phase, qRT-PCR validation studies were performed on all available exacerbation sputum samples from the whole population (n=40), and in these analyses, FEV1/FVC at enrollment, but also during the exacerbation and at convalescence, was treated as a quantitative trait.

Undetectable qRT-PCR data points were substituted with half the lowest value. Proportions were compared using Pearson chi-square or Fisher Exact test. Association was assessed using non-parametric Spearman’s Rank Order correlation coefficient. Tobit regression was used to adjust for confounders when the dependent variable was left censored. Linear regression was used to adjust for controller medication use when the dependent variable was normally distributed. T-tests, oneway ANOVA were used to assess differences between groups for normally distributed data. Mann-Whitney U, Kruskall-Wallis H were used to assess differences between groups for non-normally distributed data. Wilcoxon rank test was used for non-parametric paired data assessments. The data were analyzed using SPSS 17.0 (SPSS Inc, Chicago, IL), STATA 10.0 (StataCorp, College Station, TX) and GraphPad Prism (GraphPad Software, Inc, CA).

Results

Study population and follow-up

The study population consisted of 218 mild/moderate persistent asthmatic children aged 6–18 years who had at least one acute asthma exacerbation during the previous year. Of these 218 children, 117 experienced an exacerbation during the observation period, and sputum samples from 40 of these subjects were available for this study. Of all baseline characteristics, only ethnicity, maternal smoking, and frequency of ever having been hospitalized for asthma were different between included and excluded children who had exacerbations (Table 1).

Table 1

Characteristics of the participants who had an exacerbation and were included and not included in the present study

Characteristics	Included% (n+/n total)	Not Included% (n+/n total)	p-value
Enrollment
Ethnicity, %
Hispanic White	70.0 (28/40)	77.9 (60/77)	0.02
Non-Hispanic White	7.5 (3/40)	11.7 (9/77)
African American	2.5 (1/40)	7.8 (6/77)
Alaska Native	15.0 (6/40)	2.6 (2/77)
>1 Race	5.0 (2/40)	0.0 (0/77)
Gender, %
Female	45.0 (18/40)	44.2 (34/77)	0.3
Male	55.0 (22/40)	55.8 (43/77)
Age, yr, mean ± SD (n)	11.0 ± 3.4 (40)	10.1 ± 2.7 (77)	0.1
Height, cm, mean ± SD (n)	143.5 ± 17.5 (40)	139.4 ± 17.8 (74)	0.2
Aeroallergen skin test positive, %	73.0 (27/37)	80.3 (57/71)	0.5
FEV1/FVC Ratio, % (n)	84.8 ± 5.7 (40)	83.1 ± 6.8 (75)	0.2
Ever Eczema, %	42.5 (17/40)	53.9 (41/76)	0.3
Ever Hay Fever, %	45.0 (18/40)	52.6 (40/76)	0.6
Ever Hospitalized for Asthma, %	15.0 (6/40)	36.0 (27/75)	0.02
Medications, %
Inhaled Corticosteroids	30.0 (12/40)	36.4 (28/77)	0.5
LTRA	20.0 (8/40)	26.0 (20/77)	0.5
Advair	32.5 (13/40)	28.6 (22/77)	0.7
Parental Ever Asthma, %
Maternal	27.5 (11/40)	32.9 (25/76)	0.7
Paternal	15.0 (6/40)	15.8 (12/76)	0.9
Parental Current Smoking, %
Maternal	5.0 (2/40)	22.4 (17/76)	0.02
Paternal	35.0 (14/40)	36.8 (28/76)	0.9
Exacerbation
Time to exacerbation, days, mean ± SD (n)	114 ± 124 (40)	134.5 ± 118 (77)	0.4
Picornavirus positive, % (n)	47.2 (17/36)	43.1 (22/51)	0.8
FEV1/FVC Ratio, % (n)	83.8 ± 5.9 (36)	81.4 ± 8.8 (63)	0.2
Medications, %
Inhaled Corticosteroids	35.0 (14/40)	31.2 (24/77)	0.7
Leukotriene Receptor Antagonist	17.5 (7/40)	20.8 (16/77)	0.8
Combination therapy	37.5 (15/40)	24.7 (19/77)	0.2
Convalescent
Exacerbation to convalescent, days, mean ± SD (n)	13.6 ± 3.8 (40)	15.6 ± 7.5 (68)	0.1
FEV1/FVC Ratio, % (n)	84.8 ± 6.0 (36)	82.4 ± 6.0 (64)	0.1
Medications, %
Inhaled Corticosteroids	42.5 (17/40)	42.6 (29/68)	0.9
Leukotriene Receptor Antagonist	17.5 (7/40)	22.1 (15/68)	0.6
Combination therapy	35.0 (14/40)	33.8 (23/68)	0.9

At enrollment, modest correlations were observed between FEV1/FVC ratio and height and age (rho=−0.33, p-value =0.035 and rho=−0.37, p-value=0.020, respectively). There was no statistically significant relation between ethnicity, gender, parental asthma, current smoking, skin test reactivity, detection of picornavirus in nasal secretions, ever hospitalized for asthma, or concurrent medication use and the baseline FEV1/FVC ratio (Table S1).

Gene coexpression network analysis of exacerbation responses; variations associated with FEV1/FVC ratio

Two groups of 10 subjects each from both extreme tails of the FEV1/FVC distribution at baseline were selected (see Methods) and their mean ± SD FEV1/FVC ratios were 78.8 ± 2.1 and 91.3 ± 2.5, respectively. Total RNA was extracted from induced sputum samples obtained from these subjects during an exacerbation, and the samples were labelled and hybridized to microarrays. A holistic approach was employed to analyse the microarray data, based on reverse engineering gene network analysis.26–28 The mechanics of the network analysis algorithm are detailed in Methods; briefly, the algorithm employs a stepwise analytical process to interrogate variations in gene coexpression patterns across all the samples to identify subnets of coexpressed genes or modules. These modules contain the molecular components that execute biological functions and underpin physiological and diseased states.26–28 As illustrated in Figure 1A, network analysis of exacerbation responses in sputum identified twelve modules of coexpressed genes. To determine if any of the modules were differentially expressed in responses from subjects with or without deficits in FEV1/FVC ratios, the modules were tested with Gene Set Analysis29 – a statistical procedure that tests for the association of a set of genes with a trait of interest. This analysis provided preliminary evidence that three modules were associated with enrollment FEV1/FVC ratios (modules II, III, and IV, Figure 1B), and accordingly these modules were selected for further study.

Figure 1

Exacerbation responses are associated with baseline FEV1/FVC ratios

Gene expression was profiled by microarray in sputum samples obtained during an acute exacerbation from asthmatic children with (n=10) or without (n=10) deficits in enrollment/baseline FEV1/FVC ratios. A: Modules of coexpressed genes were identified by reverse engineering gene network analysis of the microarray data set (n=20). B: The modules were tested for differential expression in responses from subjects with or without deficits in baseline FEV1/FVC ratios. Positive or negative values on the vertical axis indicate that a module is enriched with genes that are increased or decreased respectively (based on gene-level Bayesian t.test statistics25) in responses from subjects with deficits in FEV1/FVC ratios as compared with those with normal FEV1/FVC ratios. The p-values are derived from a statistical analysis at the module-level utilizing the Gene Set Analysis test without adjustment for multiple testing.29 ** p-value < 0.01; * p-value < 0.05.

The first module (module III, Figure 1B) was decreased in subjects with low FEV1/FVC ratios as compared with those with normal FEV1/FVC ratios at baseline, and it was designated the “Th1-like/cytotoxic pathway” because it was enriched for signatures associated with Th1/IFN-γ (CCL8, CCR5, CXCL9, CXCL10, indoleamine 2,3-dioxygenase, IFN-γ, IL-12RB2, SLAM, SOCS1, STAT4, TxK)34 and cytotoxic responses (killer cell lectin-like receptor (KLRD1), IL-15, IL-21, granulysin, granzyme B, granzyme K, perforin-1). The second module (module IV, Figure 1B), “interferon signalling pathway”, was also decreased in subjects with low FEV1/FVC ratios as compared with those with normal FEV1/FVC ratios at baseline, and it was mainly composed of genes downstream of interferon signalling (interferon-induced proteins [35, 44, 44L, H1, I6, T1, T3, T5, TM3], IRF-7, Mx1, OAS2, OAS3, OASL, protein kinase R [EIF2AK2], PML, STAT2, TRAIL). The third module (module II, Figure 1B) was elevated in subjects with low FEV1/FVC ratios as compared with those with normal FEV1/FVC at baseline, and it was designated the “epithelial differentiation pathway” because it was enriched for genes which are associated with the mucociliary differentiation of airway epithelial cells (ECM1, ELF3, CLCA4, DUOX2, FAM3D, HSPB8, KRT3, KRT4, KRT16, MAL, MUC1, PPL, RCHG, SCEL, SERPINB2, SPRR1A, SPRR3).35 The Th1-like/cytotoxic, interferon signalling, and epithelial differentiation pathways contained 115, 78, and 71 genes respectively, and the complete set of genes in these pathways are listed in Table S2, Table S3, and S4. The stability of the three selected modules was tested using a leave-one-out procedure as detailed in Methods, and was greater than 80 % for the “Th1-like/cytotoxic” and “epithelial differentiation” modules. The stability of the “interferon signaling” module was lower at 63 %, however this cluster contained a highly conserved core comprising an archetypal interferon-induced gene signature (see Table S3).

Th2-related pathways were not detected in the network analyses; however, a conventional differential expression analysis of the microarray data suggested that FCER1A and IL-5 were elevated in the responses from subjects with low FEV1/FVC ratios (Figure S1). Therefore we also selected Th2-related genes (FCER1A, IL-5, IL-13) for further study.

Gene expression patterns during acute exacerbations and at 7–14 days convalescence

The microarray analyses above were based on sputum samples obtained during an acute exacerbation, thus it is not known if the gene expression signals are transient/activation-associated, or constitutive. To obtain more detailed information in this regard, we profiled expression of representative genes with known immunological functions from the identified pathways, by real time qRT-PCR in paired acute and convalescent samples, which were available for a subset of the study population (n=18). The responses for all 18 subjects combined are shown in Table 2 (see Figure 2 for illustrative examples), and these data demonstrated that representative genes in the Th1-like/cytotoxic pathway, the interferon signalling pathway, and the Th2 pathway (FCER1A, IL-5) were significantly elevated in acute versus convalescent comparisons, and there was also a trend for increased levels of IL-13 (p-value=0.051). In contrast, genes in the epithelial differentiation pathway were not significantly modulated in acute versus convalescent comparisons (Table 2). Of note, expression of the house keeping gene HMBS which was utilized to normalize the qRT-PCR data to adjust for variations in RNA quantity and quality was not different in acute versus convalescent comparisons (p-value=0.92).

Table 2

Expression of representative genes from the Th1-like/cytotoxic, interferon signaling, and Th2-related pathways are upregulated in sputum during acute exacerbations in comparison to 7–14 days convalescence.

Gene symbol	Biological function	Exa vs Conv p-valueA
Th1-like/cytotoxic pathway
CXCL10 (IP-10)	Highly expressed during rhinovirus-induced asthma exacerbations. Attracts CD4 Tcells, CD8 T cells, and iNKT cells via CXCR3.	0.0004
IFN-β	IFN-α/β induce a robust antiviral state, are essential for immunity to most viruses, andare potent inducers of cytotoxic responses.	0.0182
IFN-γ	Principal Th1 cytokine and macrophage activating factor. Essential for immunity tomycobacterium and some viruses.	0.0013
IL-2Rα	IL-2 signaling regulates immune responses via the maintenance of Treg and inductionof activation-induced cell death, promotes the proliferation and differentiation of NKand T cells, and programs antiviral CD8 T cell memory responses.	0.0009
IL-15	Role in development, growth, homeostasis, and cytotoxic function of CD8 T cells, NKand iNKT cells.	0.0006
IL-21	Enhances cytotoxic function of CD8 T cells, NK and iNKT cells. Essential for controlof chronic viral infections. Regulates IgE responses.	0.0024
GZMB (granzyme B)	Principal mediator of CD8 T cell and NK cytotoxic responses. Role in antiviraldefence.	0.0011
PRF1 (perforin-1)	Delivers cytotoxic granules containing granzyme B to induce apoptosis of target cells.Essential for immunity to some viruses.	0.0090
Interferon signaling pathway
Mx1	Induced by IFN-α, -β, -λ. Traps viral components and inhibits viral replication.	0.0007
PKR (EIF2AK2)	Induced by IFN-α, -β, -λ. Inhibits protein translation via phosphorylation of EIF2α.	0.0011
PML	Induced by IFN-α/β and IFN-γ, forms nuclear bodies and inhibits viral replication.	0.0163
CCL2 (MCP-1)	Induced by IFN-α/β. Essential for recruitment of macrophages and protection of thealveolar epithelium during respiratory viral infections.	0.0003
Epithelial differentiation pathway
ECM1	Loss-of-function mutations in ECM1 cause Urbach-Wiethe disease, a rare, autosomalrecessive disorder characterized by abnormalities in skin and mucous membranes.	0.8875
KRT4	Structural protein expressed in differentiated layers of mucosal and esophagealepithelia.	0.8192
SPRR3	Marker of squamous epithelium.	0.6243
TMPRSS11D	Trypsin-like protease expressed in bronchial epithelial cells. May play a role in airwayremodelling.	0.3550
Th2-related pathways
FCER1A	High affinity receptor subunit for IgE which is upregulated on Mo/DC during viral-driven exacerbations in children and on DC in viral-driven lung disease in mice.	0.0060
IL-5	Role in the pathogenesis of exacerbations in eosinophilic, refractory form of asthma.	0.0084
IL-13	Mediates cardinal features of asthma and viral-driven lung diseases in mice. Role inthe pathogenesis of allergen-driven late phase response in humans.	0.0512

Normalized gene expression levels as measured by qRT-PCR were compared in paired sputum samples obtained from asthmatic children (n=18) during an exacerbation (Exa) or 7–14 days convalescence (Conv). Statistical analyses by Wilcoxon signed rank test. Undetected data points were substituted for half the lowest value. See Table S5 for references.

Figure 2

Inflammatory genes are upregulated in sputum during acute exacerbations of asthma in comparison to 7–14 days convalescence

Normalized gene expression levels were measured by qRT-PCR in paired sputum samples obtained from asthmatic children (n=18) during an acute exacerbation (Exa) and 7–14 days later (Conv). Statistical analysis by Wilcoxon signed rank test. Undetected data points were substituted for half the lowest value.

Stratification of the subjects (from Table 2/Figure 2) into FEV1/FVC tertiles demonstrated that the responses for representative genes from the Th1-like/cytotoxic and interferon signaling pathways were much more consistent and/or intense in subjects with higher baseline FEV1/FVC ratios (Figure S2–Figure S6). Moreover, statistical analyses detected consistent associations between FEV1/FVC ratios at enrollment, and exacerbation responses for genes from the Th1-like/cytotoxic and interferon signaling pathways, but not from the epithelial differentiation or Th2-related pathways (Figure S7).

qRT-PCR analysis of exacerbation responses and variations in FEV1/FVC ratios

We next sought confirmation of the preliminary associations between exacerbation response patterns and enrollment FEV1/FVC ratios detected by microarray in Figure 1B, employing real time qRT-PCR analysis of all available samples from the whole population. Of note, these subjects exhibit a broad range of FEV1/FVC ratios, thus in the analyses below FEV1/FVC was treated as a quantitative trait. The non-parametric analyses illustrated in Table 3 confirmed that expression of representative genes from the Th1-like/cytotoxic and interferon signaling pathways during exacerbations was strongly and positively correlated with enrollment FEV1/FVC ratios. However, we could not demonstrate that similar correlations were observed with FEV1/FVC ratios at the time of the exacerbations, when the sputum samples themselves were obtained. There were no associations between expression of genes from the epithelial differentiation and Th2-related pathways with FEV1/FVC ratios assessed at any stage. No association was detected with FEV1/FVC ratios and expression of the house keeping gene HMBS.

Table 3

Expression levels of Th1-like/cytotoxic and interferon signaling genes in sputum during acute exacerbations are correlated with enrollment and convalescent FEV1/FVC ratios, but not with exacerbation FEV1/FVC ratios.

	FEV1/FVC Ratio
	Enrollment (n=40)		Exacerbation (n=36)		Convalescent (n=36)
Gene symbol	Rho	p-valueA	Rho	p-valueA	Rho	p-valueA
Th1-like/cytotoxic pathway
CXCL10 (IP-10)	0.44	0.004	0.10	0.6	0.48	0.003
IFN-β1	0.38	0.016	0.25	0.1	0.41	0.01
IFN-γ	0.48	0.002	0.12	0.5	0.62	<0.001
IL-2Rα	0.43	0.005	0.08	0.6	0.51	0.001
IL-15	0.45	0.003	0.04	0.8	0.53	0.001
IL-21	0.52	0.001	0.18	0.3	0.58	<0.001
GZMB (granzyme B)	0.50	0.001	0.05	0.8	0.54	0.001
PRF1 (perforin-1)	0.53	<0.001	0.22	0.2	0.52	0.001
Interferon signaling pathway
Mx1	0.47	0.002	0.07	0.7	0.55	<0.001
PKR (EIF2AK2)	0.40	0.010	0.03	0.9	0.47	0.004
PML	0.41	0.009	0.01	0.9	0.49	0.002
CCL2 (MCP-1)	0.45	0.004	0.10	0.6	0.52	0.001
Epithelial differentiation pathway
ECM1	−0.03	0.9	0.04	0.8	−0.03	0.9
KRT4	−0.15	0.4	−0.11	0.5	−0.17	0.3
SPPR3	0.05	0.8	0.10	0.6	−0.06	0.7
TMPRS11D	−0.03	0.8	0.04	0.8	−0.03	0.9
Th2-related pathways
FCER1A	−0.12	0.5	−0.18	0.3	0.15	0.4
IL-5	0.17	0.3	−0.06	0.7	0.24	0.2
IL-13	0.24	0.13	0.12	0.5	0.22	0.2
Reference gene
HMBS	−0.20	0.2	0.05	0.8	−0.09	0.6

Gene expression levels were measured by qRT-PCR in sputum samples obtained from asthmatic children during an exacerbation. Statistical analysis by Spearman Rank correlation, unadjusted. Undetected data points were substituted for half the lowest value. High quality lung function data was available for all subjects at enrollment, but was either missing or did not meet ATS guidelines for some subjects during the exacerbation or at convalescence.

Inclusion in the analyses of sputum samples which contain a higher proportion of contaminating squamous cells may decrease power to detect true associations, and this could potentially explain the lack of association between the epithelial and/or Th2-related pathways with FEV1/FVC ratios. To address this issue, we excluded samples which contained more than 30 % squamous cells, repeated the analyses, and the results were unchanged (Table S6).

Variations in gene expression could potentially be explained by variations in the clinical characteristics of the study population. However, no consistent associations were found between gene expression levels and skin test reactivity, or evidence of concurrent picornavirus infection, or use of specific medications at the time of exacerbation (Table S7).

Exacerbation responses, specific medications, and FEV1/FVC ratios

Linear regression modelling was employed to investigate associations between exacerbation response patterns and FEV1/FVC ratios, whilst adjusting for medications used at the time of the exacerbation. Consistent with the unadjusted analyses in Table 3, we identified strong and positive associations between the Th1-like/cytotoxic and interferon pathways and enrollment/convalescent FEV1/FVC ratios, and again, no associations were detected between FEV1/FVC ratios and either Th2 or epithelial differentiation pathways (Table S8).

Cellular immune signatures enriched within the Th1-like/cytotoxic pathway

In order to attempt to infer the cellular origins of the gene expression signals in the Th1-like/cytotoxic pathway, we employed a bioinformatics “in-silico” approach using a publicly available database (see Methods). Previous studies have documented high expression levels of the cytotoxic pathway in CD4 T cells, CD8 T cells, and NK cells from blood,36 and in these same cell populations isolated from the airways during respiratory viral infections.37 Our bioinformatics analyses (Table S9) showed that the list of genes in the Th1-like/cytotoxic pathway was significantly enriched for molecular signatures associated with antiviral responses (p-value = 1.13 × 10−9), CD4 T cells (p-value = 1.05 × 10−9), cytotoxic CD8 T cells (p-value = 1.63 × 10−5), T cell receptor signalling (p-value = 2.7 × 10−4), and NK cell mediated cytotoxicity (p-value = 1.77 × 10−5). However, we also identified a signature for iNKT cells (p-value = 1.04 × 10−5) in the Th1-like/cytotoxic pathway. In addition, two genes (SLAM, SAP/SH2D1A, Table S2) which are central to the development of iNKT cells38 are also part of the Th1-like/cytoxic network.

To obtain more detailed information on the role of T cells and iNKT cells in exacerbation responses, we investigated expression of T cell receptor transcripts by qRT-PCR in the sputum samples.20 Expression of the T cell receptor constant chain (TRBC2) – a marker of T cells including iNKT cells, was detected in all exacerbation and convalescent samples, and expression levels were significantly elevated during exacerbations (Figure 3A). Expression of the iNKT cell markers Vα24 and Vβ11 was detected in 80 % and 70 % of the exacerbation responses respectively (data not shown), and once again expression levels were significantly elevated during exacerbations (Figure 3B/C).

Figure 3

Markers of T cells and iNKT cells are upregulated during exacerbations; the latter are highly correlated with IL-12A and IL-21 responses

A–C: Normalized gene expression levels measured by qRT-PCR in paired sputum samples obtained from asthmatic children (n=18) during an acute exacerbation (Exa) and at 7–14 days convalescence (Conv). Statistical analysis by Wilcoxon signed rank test. D–E: Normalized gene expression levels measured by qRT-PCR in sputum obtained during an acute exacerbation from asthmatic children (n=40). Statistical analysis by Spearman Rank correlation. Undetected data points were substituted for half the lowest value.

Finally, we investigated associations between expression of T cell and iNKT cell markers, exacerbation response patterns, and FEV1/FVC ratios. These analyses demonstrated that TRBC2 and iNKT cell markers were not related to FEV1/FVC ratios (data not shown), however the former was strongly and positively correlated with the activation of Th1-like/cytotoxic and interferon pathways, and the latter were strongly correlated with IL-12A and IL-21 responses (Table S10, see Figure 3E/F for illustrative examples).

Discussion

Childhood asthma is a heterogeneous condition characterized by recurrent episodes of reversible airway obstruction. In a substantive proportion of children with asthma, the disease is associated with progressive airflow limitation, as assessed by the development of deficits in FEV1/FVC ratio, a useful index of the remodelling state of the airway unrelated to anthropometric data.22 What determines this irreversible loss of lung function is unknown, but a recent longitudinal study in children and adults with asthma of recent onset suggested that those who had exacerbations during follow-up were at increased risk.10 The purpose of this study was thus to determine if children with baseline airflow obstruction had a pattern of inflammatory responses during acute asthma exacerbations that was different from that of children without baseline airway obstruction. We demonstrated that decreased activation of Th1-like/cytotoxic and interferon signaling pathways during acute asthma exacerbations were strongly associated with chronic airways obstruction. These associations were independent of atopy, the detection of concurrent picornavirus infections, and the use of medications at the onset of the exacerbation. Th2-related pathways were also detected in the exacerbation responses (Figure 2, Figure S1), but variations in these pathways were not related to FEV1/FVC ratios. Interestingly and unexpectedly, the patterns of association with lung function observed between acute gene expression and FEV1/FVC ratios at baseline and during convalescence were not observed with FEV1/FVC ratio measured during the acute episode. We speculate that determinants of the severity of the acute obstructive response may be different from those that influence long-term effects on lung function. Finally, bioinformatics analyses of cellular immune signatures enriched within the Th1-like/cytotoxic pathway suggested a role for T cells and iNKT cells in the exacerbation responses, and we showed that the activation of Th1-like/cytotoxic and interferon signaling networks was strongly correlated with a marker of T cells (TRBC2). Markers of iNKT cells (Vα24, Vβ11) were also associated with acute exacerbations, and although they were not related to airways obstruction, they were strongly correlated with IL-12A and IL-21 responses. These results thus suggest that impairment of acute Th1-like/cytotoxic and interferon signaling responses presumably to viruses and other environmental exposures may play a major role in the development of chronic airways obstruction in asthma.

Our findings confirm and extend previous studies demonstrating deficient rhinovirus-induced interferon responses in asthmatic adults in vivo.14, 15, 17 Of particular interest were the studies by Wark et al13 who compared interferon beta responses, induction of apoptosis, and viral release after rhinovirus infection in bronchial epithelial cells obtained from asthmatic adults with and without airflow limitation (with only the former requiring treatment with inhaled corticosteroids) and normal controls. They found that, whereas interferon beta and apoptotic responses were impaired in both groups with asthma, viral release from these cells was significantly increased only among asthmatic subjects with airflow limitation. Contoli et al14 reported deficits in interferon lambda responses and increased viral shedding by rhinovirus-infected bronchial epithelial cells obtained from subjects with asthma who had mild baseline airflow limitation. These results thus suggested that deficits in several innate inflammatory responses might facilitate virus replication and cytolysis, with increased infection of neighboring airway cells, thus inducing exaggerated secondary responses, which in turn may activate remodeling and abnormal repair mechanisms. Our study addresses this potential complexity by adopting an unbiased approach to the assessment of inflammatory networks that are set in motion during real life asthma exacerbations in children. We demonstrate that the Th1-like and interferon pathways function in the broader molecular context of the cytotoxic pathway (granulysin, granzyme B, perforin-1) and common gamma chain cytokine signalling pathways (IL-2Ra, IL-15, IL-21), which are known to enhance the differentiation, proliferation, and/or cytotoxic functions of CD4 T cells, CD8 T cells, NK cells, and iNKT cells (Table 2). It is also noteworthy that interferon beta in addition to inducing a robust antiviral state via upregulation of the archetypal interferon-induced genes (Mx1, PKR, PML) also potently activates cytotoxic responses.39 Thus, activation of a complex network of interrelated immune mechanisms seems to be necessary to protect subjects with asthma from the development of airflow limitation after acute asthma exacerbations.

The detection of iNKT cells in sputum during exacerbations is not all that surprising, given their established role in immune responses against respiratory infections18, 40 and to a broad range of other pathogens.41 Our findings demonstrate that markers of iNKT cells are strongly associated with IL-12A and IL-21 responses (Table S10, Figure 3E/F). Activated iNKT cells potently induce production of the Th1-promoting cytokine IL-12 by dendritic cells,42 and IL-12 drives iNKT cell-dependent cytotoxic responses.43 IL-21 is highly expressed by CD4 T cells and iNKT cells,44 and IL-21 enhances the cytotoxic functionality of NK and iNKT cells.44, 45 In the direct context of viral infections, deficits in IL-21 responses in cord blood are associated with increased risk for severe respiratory infections.46 Moreover, studies in animal models have demonstrated that the IL-21 pathway is essential to sustain antiviral responses and control persistent infections.47 Based on the strong association of iNKT cell markers with IL-12A and IL-21 responses, we hypothesize that iNKT cells play a major role in coordinating protective Th1-like and cytotoxic responses during exacerbations.

Studies in animal models have suggested a pathogenic role for iNKT cells in asthma.18, 48 However, immune responses in mice do not always mimic those in humans,49, 50 and data on iNKT cells in human asthma are limited and conflicting.19, 20 iNKT cells were significantly increased during exacerbations in our data, but this association was modest due to evidence for the persistence of these cells in sputum 7–14 days after the exacerbation in a subset of the subjects (Figure 3B/C). This was in contrast to a broader marker for T cells, which was reduced to relatively low levels at 7–14 days convalescence (Figure 3A). Unfortunately, numbers of subjects and duration of follow up were insufficient to determine if persistence of markers for the presence of iNKT cells was associated with chronic airflow limitation, as suggested by mouse models. Further studies are needed to examine the relationship between acute asthma exacerbations, chronic airway obstruction, and numbers of iNKT cells and T cells in the airways.

This study has limitations which should be acknowledged. In our study, the abnormalities in lung function were already present at enrollment, and thus, a cause-effect relation between the inflammatory patterns we observed and these abnormalities cannot be proven using our study design. However, a study design in which these patterns were measured before the development of airflow limitation and the latter was assessed prospectively thereafter would require testing young children in whom induced sputum cannot be obtained and many years of follow-up. Also, the expression profiling studies and network analyses were performed on a mixed cell population from induced sputum, which is subject to sampling variability and contamination with saliva and squamous cells. Follow up studies on highly purified cell populations should thus increase the sensitivity and precision/resolution of these analyses. Finally, sputum induction was not performed during severe exacerbations because it was deemed unsafe and unethical, and thus the mechanisms operating in severe exacerbations may be different from those reported herein.

These limitations notwithstanding, our findings for the first time provide a global perspective of the immunological networks that are operating during naturally occurring acute asthma exacerbations in children, and broadly characterize variations in these networks which are associated with deficits in baseline FEV1/FVC ratios. Our results identify a large range of logical candidates with well characterized immunomodulatory properties (Table 2) for intervention studies to prevent the development of chronic airflow limitation in children with asthma.