Whole transcriptome analysis of high and low IFN-α producers reveals differential response patterns following rhinovirus stimulation.

OBJECTIVES: Viral respiratory infections cause considerable morbidity and economic loss. While rhinoviruses (RV) typically cause little more than the common cold, they can produce severe infections and disease exacerbations in susceptible individuals, such as those with asthma. Variations in the regulation of key antiviral cytokines, particularly type I interferon (IFN-α and IFN-β), may contribute to RV susceptibility. To understand this variability, we compared the transcriptomes of high and low type I IFN producers.
METHODS: Blood mononuclear cells from 238 individuals with or without asthma were cultured in the presence or absence of RV. Those samples demonstrating high or low RV-stimulated IFN-α production (N = 75) underwent RNA-sequencing.
RESULTS: Gene expression patterns were similar in samples from healthy participants and those with asthma. At baseline, the high IFN-α producer group showed higher expression of genes associated with plasmacytoid dendritic cells, the innate immune response and vitamin D activation, but lower expression of oxidative stress pathways than the low IFN-α producer group. After RV stimulation, the high IFN-α producer group showed higher expression of genes found in immune response biological pathways and lower expression of genes linked to developmental and catabolic processes when compared to the low IFN-α producer group.
CONCLUSIONS: These differences suggest that the high IFN-α group has a higher level of immune system readiness, resulting in a more intense and perhaps more focussed pathogen-specific immune response. These results contribute to a better understanding of the variability in type I IFN production between individuals.

Introduction

Rhinoviruses (RV) remain the most common respiratory virus causing common colds, even during the SARS‐CoV‐2 pandemic. While RV infection is seemingly innocuous in most healthy people, it can cause complications in vulnerable groups with chronic respiratory diseases such as cystic fibrosis, idiopathic pulmonary fibrosis, chronic obstructive pulmonary disease and asthma. Indeed, RV causes more asthma exacerbations than the influenza virus. RV infection worsens airway inflammation and may predispose to secondary bacterial infections, though whether these are due to abnormal host immune defence, to RV‐induced immunopathology or to both remains unclear.

Rhinoviruses‐stimulated airway structural cells and circulating leukocytes from people with asthma may produce insufficient amounts of type I interferons (IFN‐α and IFN‐β), as reviewed recently, although not all investigators have been able to confirm these findings. , IFN production appears to vary with asthma severity and between different inflammatory phenotypes, suggesting that IFN insufficiency may be characteristic of specific asthma subtypes, rather than being common to all with asthma. Recent study has linked both high and low IFN production to acute wheezing in children with respiratory illness, further supporting variable IFN responses to viruses in asthma patients.

Antiviral immune responses depend on pathogen pattern receptors, such as melanoma differentiation‐associated gene 5 (MDA‐5), retinoic acid‐inducible protein I (RIG‐1) and various Toll‐like receptors (TLRs) that recognise specific viral molecular patterns and induce type I IFN production. While all nucleated cells produce IFN‐α and IFN‐β, plasmacytoid dendritic cells (pDC) have the greatest capacity to produce these cytokines and are central to antiviral immune responses.

Optimal IFN‐α and IFN‐β production is critical for the antiviral immune response, as too little or too much IFN has adverse health consequences. Hence, type I IFN production needs to be tightly regulated, involving multiple checkpoints and complex patterns of gene expressions. At a molecular level, the interferon regulatory factor (IRF)7 is a major hub regulating type I IFN production, while at a cellular level, pDC are recognised for their capacity to rapidly produce large amounts of IFN‐α, with recent publications highlighting key mechanisms. ,

Although such studies have shed important light on type I IFN biology, it is also important to understand how IFN‐α production varies between individuals. Investigating differences between people exhibiting efficient and weak antiviral responses may provide important insights into infection susceptibility. Differences occurring before and during infection may contribute to efficient/weak antiviral responses. While a variety of researchers have made important observations concerning antiviral immune response variability, , , , , , , these studies have not examined RV, which is known to induce distinct host responses relative to other viruses. To address this important knowledge gap, the primary study aims were to define the folllowing:

In the absence of RV stimulation which transcriptional variations differentiate high IFN‐α and low IFN‐α producers, and

During RV stimulation which transcriptional variations differentiate high IFN‐α and low IFN‐α producers.

The secondary study aim was to determine whether transcriptional variations differed between those with asthma and healthy participants in unstimulated and RV‐stimulated conditions.

We made use of samples obtained during our recent study of the immunological and clinical variables associated with cold frequency in 301 individuals. RV‐stimulated peripheral blood mononuclear cells (PBMC) showed a very broad range of IFN‐α production in these individuals. From these samples, we selected a subsample of 75 individuals with contrasting high IFN‐α and low IFN‐α production in order to perform differential gene expression analysis in these ‘extreme phenotypes’. Because the study groups differed in mean body mass index (BMI), statistical methods were used to assess whether BMI might be affecting gene expression.

Results

Marked variability in RV‐stimulated IFN‐α production

Cells cultured without RV stimulation elicited no detectable IFN‐α production (below the 25 pg mL−1 detection limit), whereas RV stimulation of PBMC for 24 h elicited a very broad range of IFN‐α production across 238 participants with European ancestry (Figure 1; median 1005.2 pg mL−1; IQR 560.7, 1690.6). In order to provide more detailed analysis of gene expression patterns, we focussed on the extreme phenotypes of IFN‐α production: 75 participants with the highest IFN‐α and lowest IFN‐α production were selected from asthma and healthy groups to contrast differences in whole‐genome transcriptomes (Figure 2). The analysis was restricted to participants with European ancestry, in order to minimise confounding by ancestral background. Four samples did not pass quality control during RNA sequencing or principal component analysis (PCA) including one control asthma IFN‐α‐low sample, one control healthy IFN‐α‐high sample, one RV healthy IFN‐α‐high sample and one RV healthy IFN‐α‐low sample. These four samples were excluded from RNA sequencing, thereby reducing the total number of samples used in the complete analysis from 150 to 146 (Supplementary table 1). Further details regarding study design are shown in Figure 2.

Figure 1

Variation in RV‐stimulated IFN‐α production. Frequency distribution of RV‐stimulated IFN‐α production in 238 PBMC samples tested in two biological replicates. RV‐stimulated samples with low IFN‐α production (lowest 15% of IFN‐α production) produced less than 414 pg mL−1 IFN‐α (orange‐shaded), and samples with high IFN‐α production (highest 15% of IFN‐α production) produced more than 2200 pg mL−1 IFN‐α (green‐shaded). The black dotted line presents the frequency distribution fitted to a Gaussian distribution. The y‐axis shows relative frequency of counts in each IFN‐α value range.

Figure 2

Study design. Of the 301 study participants, 238 had self‐identified European ancestry. Samples for gene expression analysis comprised those in the lowest 15% and highest 15% IFNα production within the asthma group and those within the lowest 15% and highest 15% IFNα production within the healthy group. When RNA quality was suboptimal, samples did not proceed to RNAseq.

Clinical and molecular features of the high and low IFN‐α participants

The only criteria used to select samples for transcriptome analysis were high or low IFN‐α production (Figure 2) and the presence or absence of asthma. It is noteworthy that the IFN‐α‐low groups had higher BMI (24.0 and 27.0 kg m−2 vs 23.0 and 24.0 kg m−2; P = 0.008) than the IFN‐α‐high group (see Supplementary table 1). There were age differences between the four groups (P = 0.042): the asthma IFN‐α‐low group was oldest (median age 47 years), whereas the asthma IFN‐α‐high group was youngest (median age 27 years) with the two healthy groups in between (healthy IFN‐α‐high median 32 years; healthy IFN‐α‐low median 34 years). Minor differences in sex ratios were observed between groups, but these were not statistically significant. There were no significant differences in asthma control, asthma severity or self‐reported respiratory infection frequency between IFN‐α‐high and IFN‐α‐low groups (Supplementary table 2).

Gene expression signatures of IFN‐α producer groups are not determined by BMI

As the IFN‐α‐low groups had significantly higher BMI than the high IFN‐α producer groups, this raised a concern that BMI differences could confound the analysis. We used statistical two‐way ANOVA testing to identify which differentially expressed genes (DEGs) between high and low IFN‐α producer groups were likely confounded by BMI. The participants were divided into two groups with BMI below or at/above 25 kg m−2, the threshold for normal weight and overweight condition. The resulting interaction between IFN‐α producer groups and BMI groups identified the ability of BMI to confound the signatures of IFN‐α producer groups. We observed that only a small fraction of genes show a significant interaction between BMI and IFN‐α producer groups (baseline – 2%, RV treatment – 5%) and are likely to behave differently in obese and normal weight samples; however, majority of genes are used reliably to differentiate gene expression between high and low IFN producer groups (Supplementary figure 1). Further details on the DEGs that associate with BMI regardless of IFN‐α producing ability are listed in Supplementary figure 2.

High IFN‐α producers have more circulating plasmacytoid dendritic cells and produce more TNF

As pDC are the most potent IFN‐α‐producing cells in the circulation, variations in blood pDC numbers are likely to impact on IFN‐α production in response to RV stimulation. Hence, we compared pDC quantity and cytokine production by RV‐stimulated PBMC across the 75 control samples stratified by IFN‐α production and presence or absence of asthma. We have previously shown that the gene expression of C‐type lectin domain family 4 member C (CLEC4C) can be used to quantify circulating pDC in whole blood. As expected, CLEC4C gene expression in whole blood was significantly lower in the IFN‐α‐low groups signifying lower quantities of the IFN‐α‐producing pDC (Figure 3; P‐value < 0.001).

Figure 3

CLEC4C gene expression in whole blood (a) and cytokine production in PBMC by study group (b, c). Individual datapoints are overlaid with median (bold line) and interquartile range markers (thin lines) for each group. Statistically significant Kruskal–Wallis ranks sum tests between the four groups and significant post hoc Mann–Whitney U‐tests are indicated as ***P < 0.001; **P < 0.01; *P < 0.05.

We also assessed production of other cytokines by RV‐stimulated PBMC. The production of the TLR8‐stimulated pro‐inflammatory cytokine TNF was also lower in the IFN‐α‐low groups (Figure 3; P‐value = 0.006), whereas the production of another TLR8‐activated pro‐inflammatory cytokine IL‐12 was similar in all four groups.

Whole‐genome transcriptomes in different study groups

The whole‐genome transcriptomes were compared across the 75 participants stratified by IFN‐α production and the presence or absence of asthma. PCA transcriptomes indicated a marked difference between unstimulated and RV‐stimulated samples (not shown), and a modest grouping of samples by IFN‐α producer groups across PC1 [Figure 4 unstimulated (a), RV‐stimulated (b)]. Because PCA revealed no apparent distinction between asthmatic and non‐asthmatic samples, it seemed reasonable to undertake a pooled analysis of the IFN‐α‐high and IFN‐α‐low groups, regardless of asthma status, to address factors associated with variations in IFN‐α production – the primary objective of the study.

Figure 4

PCA of samples. (a) Unstimulated samples and (b) RV‐stimulated. Green indicates high and orange indicates low IFN‐α producer samples. Round shapes indicate asthma and triangle shapes healthy samples. Three samples (one unstimulated and two RV‐stimulated) were found suboptimal and were excluded from subsequent analyses.

Baseline gene expression profiles differ between high and low IFN‐α producer groups

We next performed a separate analysis of the transcriptomes using unstimulated and RV‐stimulated PBMC from the IFN‐α‐high and IFN‐α‐low producer groups. With a cut‐off of false discovery rate (FDR) q‐value < 0.05 and log fold‐change > 1 of gene expression, unpaired differential gene expression analysis in the unstimulated samples returned 39 genes with high expression (Supplementary table 3) and 10 genes with low expression (Supplementary table 4) in the IFN‐α‐high group compared to the IFN‐α‐low group. The heatmap visualisation in Figure 5 reveals a major subcluster of genes with high expression that have consistently high expression in the IFN‐α‐high group and low expression in the IFN‐α‐low group.

Figure 5

Differential gene expression in 73 unstimulated samples. Gene expression of the DEGs in unstimulated PBMC samples IFN‐α‐high group vs IFN‐α‐low group. Log2‐transformed gene expression is presented as a difference from the median red denoting high expression and blue low expression. Samples are hierarchically clustered with the Ward method and genes with the weighted pair group method with arithmetic mean. Samples are colour‐coded for asthma status and IFN‐α producer group.

Sample clustering with Ward hierarchical cluster analysis positions asthma samples randomly, which supports our results from the PCA (Figure 4) that asthma status does not have a major impact on differential gene expression in our samples. We searched the manually annotated records in the protein knowledgebase UniProt and PubMed (NIH, USA) for any additional literature supporting the understanding of the functions of the DEG products. Common functional groups that arose from the literature are annotated for each gene in Supplementary tables 3 and 4.

In comparison with the IFN‐α‐low group, the genes with high expression in the IFN‐α‐high group were most often related to immune function, whereas the genes with low expression were most often related to oxygen transport. The immune‐related genes included those associated with the complement system (C1QA, C1QB, C1QC, C3, CD93), antigen presentation (CLEC4C, CPVL, RNASE1, ASGR2), B‐cell chemokine (CXCL13), immune regulation (IDO1, CHI3L1, UBD, VCAN), macrophage and pDC surface markers (CLEC4C, MS4A4A, MS4A6A). Three potassium channel encoding genes (KCNJ10, KCNJ15 and KCNMA1) and a transcription factor gene ETV3L associated with vitamin D activation also had high expression in the IFN‐α‐high samples.

Several genes with high expression in IFN‐α‐high have been associated with allergy [CD93 ; AOC1 ; FXYD6 ], asthma risk [CHI3L1 ], airway remodelling [VCAN ; KRT5 ] and airway obstruction [PALD1 ]; however, no distinct pattern in asthma samples can be shown in relation to these genes.

A group of oxygen transport genes that have low expression encodes haemoglobin subunits (HBA1, HBA2, HBB, HBG2) and an enzyme involved in haem biosynthesis (ALAS2). The gene group forms a subcluster with low expression in most IFN‐α‐high samples. Four other genes with low expression have a role in wound healing (MET, DPYSL3), stress response (HSPA1B) and as inflammatory markers (HSPA1B and NR4A2).

RV‐stimulated gene expression differs between high and low IFN‐α producers

We next assessed differential gene expression in RV‐stimulated PBMC. With a cut‐off of FDR q‐value < 0.05 and log fold‐change > 1, unpaired differential gene expression analysis returned 55 genes with high expression (Supplementary table 5), and 73 genes with low expression (Supplementary table 6) in the IFN‐α‐high group. Like the clustering analysis of the unstimulated samples (Figure 5), the RV‐stimulated PBMC samples show expression profiles that cluster primarily by IFN‐α‐producer group, rather than by the presence or absence of asthma as presented in Figure 6.

Figure 6

RV‐stimulated differential gene expression. Gene expression of the DEGs in 73 RV‐stimulated PBMC samples IFN‐α‐high group vs IFN‐α‐low group. Log2‐transformed gene expression is presented as a difference from the median red denoting high expression and blue low expression. Samples are hierarchically clustered with the Ward method and genes with the weighted pair group method with arithmetic mean. Samples are colour‐coded for asthma status and IFN‐α producer group.

This clustering analysis reveals three distinct clusters of genes with similar expression patterns (Figure 6). The gene cluster in the middle contains the upregulated genes, including IFN genes. The expression of those genes is uniformly high in the IFN‐α‐high group and low in the IFN‐α‐low group.

In contrast, the genes with low expression separate into two main clusters at the top and bottom of the heatmap. The expression of those genes is consistently low in the IFN‐α‐high group samples but is more varied in the IFN‐α‐low group samples. The top gene cluster contains several subclusters of genes with high expression in distinct IFN‐α‐low sample clusters. In contrast, the entire bottom cluster of genes has high expression in only one IFN‐α‐low sample cluster consisting of mostly non‐asthma samples. Some clustering with asthma and healthy samples is evident, and these sample subclusters have a specific gene subcluster expression pattern.

We searched the manually annotated records in the protein knowledge database UniProt and any additional literature available for the function of the DEG products. The most common relevant functions are annotated for each gene in Supplementary tables 5 and 6.

As expected, IFNB1, IFNL1 and IFNW1 had high expression in the IFN‐α‐high samples (Supplementary table 4). Including those genes, 33 out of the 56 most genes with high expression are related to immune function, and two have supporting roles in the antiviral immune response: breakdown of viral lipid membrane [ELOVL7 ] and nucleic acid metabolism (UPB1). Interestingly, several genes with low expression in IFN‐α‐high group have antibacterial functions, including BPI, HP, HAMP, LTF, LCN2, PGLYRP1 and STAB1.

The other immune‐related genes with high expression have a known cytokine (IL18, IL31RA, CXCL9, FLT), antibacterial [DEFB1; ACOD1; AQP9 ], antiviral (DEFB1), chemokine (CCL18, CCL19, CCL23, CXCL13) or antigen presentation (CD1D, FCGR1B) function.

Conversely, under RV stimulation, cells from the IFN‐α‐high group showed low expression of gene products related to structural or extracellular matrix (ECM) function. Three matrix metalloproteinases (MMP9, MMP7, MMP8) involved in ECM breakdown were prominent, while the matrix metalloproteinase inactivator TIMP3 showed only low expression. SERPINE1 and COL23A1 are also associated with ECM organisation. The gene products of COL23A1, GREM1, STAB1, CEACAM8, CEACAM6, DSC1, FLRT2, ITGA11 and TIMP3 have a role in cell‐cell or cell‐ECM adhesion. Similarly, SEMA4C, PDGFC, MMP9, MMP7, MMP8 and CEACAM6 upregulate cell migration, whereas PODN, CYP1B1, SERPINE1 and DPYSL3 downregulate cell migration. The gene products of PDGFC, CEACAM6 and GREM1 upregulate cell proliferation, whereas PODN and CD9 downregulate cell proliferation. Also, GPC4, MET and SDC2 may be involved in cell proliferation.

Three genes had high expression in both unstimulated and RV‐stimulated samples, while seven genes showed only low expression. UBD, which was highly expressed, encodes an ubiquitin‐like protein that tags proteins for degradation and regulates TNF‐mediated NFKB signalling and dendritic cell (DC) activation/maturation. CXCL13 and PPFIA4 genes were highly expressed; they encode a B‐cell chemokine and focal adhesion regulator, respectively. IFN‐α‐high samples showed low expression of the haemoglobin subunits (HBA1, HBA2, HBB, HBG2) and haem biosynthesis enzyme (ALAS2), cytoskeleton remodelling (DPYSL3) genes and MET, a gene encoding a transmembrane receptor involved in cellular proliferation and fibrosis.

Biological pathway analysis comparing differentially expressed genes in the IFN‐α high and low producers

Gene set enrichment analysis (GSEA) was used to search for functional enrichment within the DEG sets. Figure 7a presents the gene ontology (GO) pathways associated with DEGs from unstimulated samples for the IFN‐α‐high group compared to the IFN‐α‐low group. The pathways with high expression predominantly reflect upregulation of multiple immune functions, whereas the pathways with low expression are much more varied in their function. A network representation visualises shared genes in selected pathways (Figure 7c), demonstrating considerable gene overlap within both up‐ and downregulated pathways.

Figure 7

Gene ontology pathways. (a, b) Gene ontology pathways that are enriched in the set of DEGs in unstimulated (a) and RV‐stimulated (b) PBMC samples of IFN‐α‐high group in contrast with IFN‐α‐low group. FDR, false discovery rate; NES, normalised effect size; size, number of genes found in the pathway. Pathways are sorted by significance on the horizontal axis and effect on the vertical axis as indicated by the dots. Size of the dot is relevant to the number of genes found in the pathway and colour to the effect size. (c, d) Associations between selected downregulated (blue) and upregulated (red) gene ontology pathways in unstimulated (c) and RV‐stimulated (d) PBMC samples. Pathways are enriched in the DEGs in the samples of IFN‐α‐high group in contrast with IFN‐α‐low group. Pathways with least overlap in function were selected. Line thickness reflects the number of shared genes.

Several downregulated pathways were enriched in the list of DEGs. Prominent among these important pathways were oxidative stress response genes, with high statistical significance, driven by a relatively small number of genes, including haemoglobin subunits (HBA1, HBA2, HBB) and inflammatory markers (HSPA1B and NR4A2). Together, the downregulated gene functions and the enriched GO pathways indicate a differential state of cellular stress between the IFN‐α producer groups, such that the cells from high IFN‐α producing individuals showed efficient downregulation of oxidative stress gene expression, whereas a subset of samples from low IFN‐α producing individuals exhibited high expression of oxidative stress genes. The oxidative stress pathway was not enriched in the RV‐stimulated samples.

As expected, the significantly upregulated GO pathways during RV stimulation relate to activation of antiviral immune responses (Figure 7b), whereas the genes with low expression are involved in cell proliferation, cell‐cell communication and cell migration, which is reflected in the enriched pathways.

A selection of a representative up‐ and downregulated pathways is presented in a network shown as Figure 7d. The RV‐stimulated pathways show high interconnectivity and several shared genes between the upregulated pathways, whereas the downregulated pathways share less genes and are less connected.

Validation of the DEGs using microarray results from a separate cohort

To validate the results from the differential gene expression analysis, we accessed a similar cohort of RV‐stimulated PBMC samples from a separate cohort of 17 asthma samples and 17 controls collected in 2014, as previously reported. The comparison between the validation dataset and the current dataset is shown in Supplementary tables 1 and 2.

We used the microarray results from the validation cohort to confirm the association between the expression of the DEGs and IFN‐α production. Linear regression was used to test associations between gene expression and IFN‐α production as IFN‐α production in the validation dataset was lower than in the current study, and there were insufficient numbers to restrict the analysis to high and low IFN‐α producers. Because of lower numbers of genes analysed by the microarray experiment, only 11 out of the 45 unstimulated DEGs from the main study were detected in the validation microarray (Figure 8a), and 36 out of the 129 RV‐stimulated DEGs were detected in the validation microarray (Figure 8b). After correcting for multiple testing, the unstimulated DEG CLEC4C gene expression was validated as significantly associated with IFN‐α production in the validation cohort and RV‐stimulated DEGs IFNW1 and GCOM1 gene expression was validated as significantly associated with IFN‐α production (Supplementary table 7).

Figure 8

DEGs from the main study that are detected in the validation dataset. Unstimulated condition diagram (a) and RV‐stimulated condition diagrams (b) show the DEGs that are significantly associated with IFN‐α production in the validation cohort.

Discussion

The primary aim of this study was to examine transcriptional variations associated with RV‐induced IFN‐α production. Having shown the broad range of IFN‐α production in a large group of study participants, we focussed on the extreme phenotypes, contrasting whole‐genome transcriptomes in the 15% with the highest or the 15% with the lowest IFN‐α production. While others have comprehensively profiled variability in antimicrobial immune responses, a strength of the current study was the large cohort in which weak and strong IFN‐α producers were identified prior to transcriptomic analysis. Previous studies have not comprehensively examined variations in the RV‐induced immune response, which is particularly important, given the frequency with which RV induces exacerbations of asthma and other respiratory illnesses and its predominance relative to other common respiratory viruses.

Contrary to our expectation, we saw only modest differences in gene transcription patterns between those with asthma and healthy participants. A subgroup of people with both asthma and low IFN‐α production showed high expression of two gene clusters (Figure 6). This supports the notion that only a subgroup of people with asthma have specific alterations in IFN‐α production and antiviral immunity, rather than this being a general characteristic of all people with asthma. ,

The unstimulated samples were used as a baseline measurement with the potential to reveal biological factors predicting the strength of IFN‐α response. The high and low IFN‐α‐producing groups demonstrated multiple DEGs, enriched in immune‐related genes and GO pathways. The DEGs expressed by specific immune cells may point towards differences in the function of those cells or merely quantities. Conventional DC are essential for detecting invading pathogens and activating and directing the adaptive immune response, while pDC have a specialised role as potent IFN‐α producers. , The pDC‐associated gene CLEC4C regulates antigen presentation and type I IFN production, while its consistent expression profile enables its use for pDC quantification. CLEC4C had higher expression in unstimulated samples in the IFN‐α‐high group than in the IFN‐α‐low group in whole blood samples and in the PBMC samples used for whole transcriptome analysis (Figures 3 and 5). This was confirmed in the validation study in which CLEC4C was the only DEG in unstimulated samples significantly associated with IFN‐α production (Figure 8a). Thus, CLEC4C expression appears as a key variable predicting the magnitude of subsequent RV‐induced IFN‐α production.

Other antigen‐presenting cell (APC)‐associated genes, including CD93, CPVL, IDO1, MS4A4A and VCAN, also had higher expression in the high IFN‐α producer group in unstimulated samples supporting the close relationship between APC and an efficient IFN‐α response. While MS4A4A is found on the surface of alternatively activated M2 macrophages, MS4A4A also forms part of a type I IFN signature in early rheumatoid arthritis. Thus, our finding of high MS4A4A expression in the high IFN‐α producer group may be a consequence of systemic inflammatory signals, rather than fitting into a simple M1/M2 paradigm. The endosome protein encoded by ASGR2 is important for APC function for sampling the extracellular environment by mediating the endocytosis and lysosomal degradation of glycoproteins. Moreover, the complement system is increasingly recognised to play an important role in antiviral immunity, and many of the genes with high expression in the IFN‐α‐high group were members of the complement system whose genes are enriched in the innate immune response GO‐pathway. Notably, others have reported that this GO pathway contains multiple genes whose expression shows a high degree of inter‐individual variability.

The RV‐stimulated samples were used to examine RV infection in vitro. A triggered antiviral response results in a substantial increase in the expression of genes involved in defence against the pathogen. The larger amount of type I IFN produced in response to RV and the IFN genes with higher expression in the IFN‐α‐high group compared to IFN‐α‐low group reflect that. The higher expression of multiple Th1 and antiviral immune response genes in the IFN‐α‐high group than the IFN‐α‐low group further indicates a more extensive antiviral response. Notably, the IFN‐α‐high group also showed higher TLR8‐induced TNF production. Given that TLR8‐induced cytokine production is dominated by TNF and IL‐12 production, rather than IFN‐α production, this suggests the IFN‐α‐high group has a greater capacity to respond to viruses and viral nucleic acids that extends beyond IFN‐α production. The difference in developmental and metabolic process pathways between the two groups also supports the notion the IFN‐α‐high group devotes their transcriptome to a robust antiviral immune response and suppresses non‐critical pathways, whereas the IFN‐α‐low group transcriptome retains the transcription of those other functions.

Of the 56 genes with higher expression in the IFN‐α‐high group, ten IFN and three non‐IFN genes correspond to significant cytokine expression quantitative trait loci (eQTL) genes. One of these genes with high expression was IL18, a TLR8‐induced cytokine that activates type 1 cytokine production in natural killer cells. The association of IL18 with immune variation has previously been documented: an IL18 eQTL modulates influenza‐induced IFNβ production by DC ; IL‐18 binding protein (IL‐18BP) restricts IL‐18 availability and is inversely associated with cytokine production. During experimental RV infection, IL‐18 appears protective, with IL‐18 concentrations in nasal and bronchial lining fluid inversely proportional to infection severity, both in healthy people and those with asthma. Our own findings described herein support the conclusion that IL‐18 production is an important component of effective host defence against severe rhinovirus infections and offers an attractive target for future research.

The high IFN‐α producer group showed lower RV‐associated expression of multiple antibacterial genes such as the antibacterial peptides bactericidal permeability increasing protein (BPI), hepcidin antimicrobial peptide (HAMP) and lactotransferrin (LTF). The IFN‐α‐low group upregulated genes for lactotransferrin and several haem biosynthesis and haemoglobin components, which are important in binding iron. Bacteria require iron for survival, and binding iron to those molecules is a powerful antibacterial mechanism. We speculate that this might reflect an immune response that is unspecific to pathogen type, whereas high IFN‐α producers respond with a more highly focussed antiviral response. There are clearly complex interactions between viruses and bacteria in the lung. Rhinoviruses are known to induce degradation of antimicrobial peptides and impair macrophage antibacterial responses. ,

Outside the immediate immunological functions, cells from high IFN‐α producers also showed efficient downregulation of the oxidative stress pathway at baseline, whereas a subset of low IFN‐α producers exhibited upregulation of the oxidative stress pathway. Oxidative stress typically results from external factors such as cigarette smoke or intrinsic factors including reactive oxygen species (ROS) produced in the context of inflammation. Exogenous oxidants can reduce IFN‐α production by pDC, while mitochondrial reactive oxygen species can inhibit TLR7 function , ; hence, oxidative stress might be a factor constraining virus‐induced IFN production. Oxidative stress is known to be present in asthma and correlates with clinical severity, , which might contribute to the weak RV‐induced IFN‐α production reported in some studies. Addressing oxidative stress in asthma might be beneficial by promoting a stronger antiviral immune response.

The groups showed differential expression of genes related to nutrient availability and function. The high IFN‐α producer group showed higher baseline expression of ETV3L, a vitamin D‐associated transcription factor and higher RV‐stimulated expression of CYP27B1, and the protein encoded by this gene converts 25‐OH vitamin D3 to its active form. The expression of these genes indicates that the vitamin D availability may be necessary for a strong antiviral response. In contrast, TCN1, a vitamin B12 binding protein, was downregulated in IFN‐α‐high RV‐stimulated samples. By downregulating the B12 binding protein, high IFN‐α producers may be making more B12 available.

Interestingly, Khoo and colleagues identified molecular phenotypes with some similarities to those identified herein in their study of upper airway specimens collected from children with acute asthma exacerbations. In their cluster analysis, they identified two distinct molecular phenotypes, one characterised by high IRF7 and high IFN expression, and a second characterised by low IRF7 expression, growth factor signalling and downregulation of IFN. Clinical features differed between these two phenotypes. There are a number of differences in study design, and Khoo and colleagues measured gene expression in upper airway cells collected from acutely unwell children without in vitro stimulation, whereas our study assessed circulating immune cells collected when participants were relatively well; gene expression was assessed in both stimulated and unstimulated cells. Nonetheless, we think that the two studies provide important complementary findings.

Association tests between gene expression and IFN‐α production in the validation dataset confirmed CLEC4C as an important baseline gene and confirmed that RV‐stimulated IFNW1 and GCOM1 are significantly associated with IFN‐α production. GCOM1 codes for GRINL1A complex locus 1 that interacts with the N‐methyl D‐aspartate (NMDA) receptor in the nervous system. There is some evidence that NMDA‐type glutamate receptors are expressed on lymphocytes and neutrophils and that its activation has functional consequences. In the central nervous system, interactions have been described between type I IFN function and NMDA receptors, though whether this is relevant in relation to RV requires further investigation. The GCOM1 gene has been reported to be involved with transcription elongation, and because of its nature as a complex transcription unit, several transcriptional variants are produced that could potentially have a wide range of functions.

The study groups were not matched for BMI, age and gender, which may have influenced the findings in this study. However, our analysis indicated that only a handful of genes showed a significant interaction between IFN‐α producer group and BMI. Piasecka et al. showed that age influences influenza‐induced IFN‐α production, but we found little evidence in the current study that age was affecting RV‐induced IFN‐α production, although we acknowledge that the participants did not cover a wide age range. Sex hormones influence antiviral immunity , and we recently observed that age is associated with respiratory infection frequency in women but not men, implying a role for sex hormones in antiviral immunity. Since the ratio of men and women in each of the groups was similar in that study and the current study, we think it is unlikely that gender imbalance had a major influence on our key findings.

For the differential gene expression analysis, we chose to combine the asthma and healthy groups because IFN‐α production was similar in both groups and asthma was not a significant factor in our principal component analysis. Despite this, our gene expression findings are in line with previous studies of healthy individuals. , Finally, to validate our results, we accessed a similar cohort of asthma and control samples gene expression data performed with microarray technology. However, gene matches to DEGs were limited because of the restricted sensitivity of the microarray technology compared to RNAseq. Because of the study design, it would have been impractical to collect a second cohort with high and low producers, and as the validation samples are closer to the IFN‐α‐low group in their IFN‐α production than IFN‐α‐high, the number of false negatives may have increased in the validation results.

Collectively, these results show that at baseline, high IFN‐α responders express essential components of the innate immunity (especially pDC‐related CLEC4C) at a higher level than the low IFN‐α responders. During the antiviral immune response, high IFN‐α producers devote more of their transcriptome to the production of antiviral cytokines and effector proteins while downregulating expression of genes encoding antibacterial proteins and cellular processes, such as development and metabolism pathways. Differences between the high and low IFN‐α antiviral response can be partly attributed to host factors regulating oxidative stress and availability of vitamin D and B12. This suggests that the high and low IFN‐α groups have varying levels of immune system readiness and capacity to translate this into a pathogen‐specific immune response. Our findings described herein contribute to a better understanding of the inter‐individual variability in type I IFN production in the context of respiratory virus infections, vaccination, asthma and autoimmune disease.

Methods

Participants

An asthma case–control study was designed to recruit 300 participants, but actually recruited 301. Details of the study cohort and participant characteristics were described recently. The experiments described herein relate to PBMC samples from the subset of 238 participants with European ancestry (Figure 1), selected in order to minimise genetic variations related to ancestry. RV‐induced IFN‐α production was measured and two contrasting IFN‐α producing groups were selected for more detailed analysis of gene expression, comprising the samples from the highest 15% or lowest 15% of IFN‐α production in the European subset. Samples in the high IFN‐α producer group were obtained from 18 participants with asthma and 19 healthy participants, while samples in the low IFN‐α producer group were obtained from 19 participants with asthma and 19 healthy participants, comprising a total of 75 participants. The high and low IFN‐α producer groups will henceforth be referred to as IFN‐α‐high and IFN‐α‐low groups.

To validate the results, we used samples from a separate cohort recruited in 2015. We had access to samples from 18 healthy participants and 17 participants with asthma from this latter study. One healthy sample that was part of both cohorts was excluded for the validation test reducing the number of healthy samples to 17. Characteristics of the two cohorts are presented in Supplementary table 1. Both studies received ethical clearance from the University of Queensland (project 2008000037) and Metro South Human Research Ethics Committees (Reference HREC/07/QPAH/146), and all participants gave written informed consent. The raw data are available from Gene Expression Omnibus repository (accession number GSE99858).

Cell cultures and ELISA

Peripheral blood mononuclear cells were isolated and freshly stimulated with RV16 or the TLR8 agonist VTX‐2337 (1 μm; Sapphire Bioscience, Waterloo Australia) for 24 h as described in detail elsewhere. , Control cells were cultured in media alone with no added stimuli. IFN‐α, TNF and IL‐12 concentrations in culture supernatants were measured with enzyme‐linked immunosorbent assay (ELISA; pan‐specific IFN‐α, Mabtech Ab, Sweden; IL‐12 (p70), BD OptEIA, BD Biosciences, USA; and TNF BD Biosciences, USA) as described recently.

Gene expression quantification

For the current study, RNA in PBMC samples was preserved in RNAprotect™ (Qiagen, Hilden, Germany) and stored at −80°C. Samples were extracted for total RNA and sequenced at Macrogen Inc. (Seoul, South Korea) with TruSeq mRNA kit using the NovaSeq6000 platform, with a minimum of 20 million, 100 bp paired‐end reads per sample. Two PBMC RNA samples were sequenced per individual: one unstimulated (baseline) and one RV‐stimulated. The data are available from Gene Expression Omnibus repository (accession number GSE166292).

The Genome Informatics Group at QIMR Berghofer, Queensland, Australia, aligned the sequence reads. Sequence reads were trimmed for adapter sequences using Cutadapt [version 1.11 ] and aligned using STAR [version 2.5.2a ] to the GRCh37 assembly with the gene, transcript and exon features of Ensembl [release 89 ] gene model. Quality control metrics were computed using RNA‐SeQC [version 1.1.8 ], and gene expression levels were quantified using RSEM [version 1.2.30 ].

Gene expression of the pDC marker CLEC4C was quantified in whole blood for the main cohort using relative quantitative RT‐PCR as described.

RNA from the validation samples was extracted, and transcriptional profiling was performed using Illumina Human HT‐12 microarrays (SanDiego, CA, USA) at the microarray facility at the Diamantina Institute, University of Queensland as described. Both RV and unstimulated PBMC samples were profiled. Illumina BeadStudio summary probe and summary control probe profiles were read into R using the read.ilmn() function available in the limma package. Background correction and normalisation were performed using the neqc() function available in limma, which uses the normal exponential convolution model for background correction followed by quantile normalisation. Distribution of probes for each array was assessed using boxplots, before and after normalisation, and no outlying arrays were identified. Probes that are not expressed were filtered out prior to analysis using a threshold of probes that are expressed in at least three arrays according to the detection of P‐values of 5%. 44 samples and 10 907 probes were included for further association analysis.

Differential gene expression analysis

EdgeR package was utilised for the discovery of DEGs in the IFN‐α‐high group compared with the IFN‐α‐low group. Library size was corrected using counts per million (CPM), which involves dividing each sample gene count by the total number of mapped reads. Trimmed mean of M‐values was used to normalise differences in RNA composition between samples with the function calcNormFactors() from the edgeR package. Function glmQLFit() fits a quasi‐likelihood negative binomial generalised log‐linear model to count data to identify DEGs. Significant genes were filtered with multiple testing corrected FDR < 0.05 and log fold‐change > 1. RV‐stimulated and unstimulated samples were tested separately.

Biological pathway analysis

Identification of significantly enriched biological pathways was performed using GSEA focussing on DEGs between IFN‐α‐high and IFN‐α‐low groups in unstimulated and RV‐stimulated samples. Fold change in gene expression was used for ranked list of input to GSEA, and GO pathways with FDR < 0.25 were considered as biological pathways significantly related to the DEGs as recommended by GSEA.

Statistical analysis

All statistical analyses were performed with R (version 3.4.4 ). Variables were tested for normality and consequently treated as nonparametric variables. CreateTableOne() function in the R package tableone was used to create the demographics table of the study group and perform statistical tests. The difference in sample distribution was tested with a nonparametric Mann–Whitney U‐test. Linear regression was used to test association. A P‐value < 0.05 was considered statistically significant after multiple testing adjustment where necessary. To determine the genes most likely associated with differences in BMI, two‐way ANOVA was first used to exclude interactions between IFN‐α producer groups and groups with BMI below/above 25 kg m−2. Then, least absolute shrinkage and selection operator (LASSO) regression analysis was used to test BMI association with gene expression.

Graphs

Boxplots and before–after plots were created with GraphPad Prism software (version 7, La Jolla, USA). GO‐pathway networks were generated in Cytoscape interface [San Diego, USA ], and other graphs were created with R using packages ggfortify, ggplot2, Pheatmap and dendsort.

Conflict of interest

JW Upham received an unrestricted grant from AstraZeneca that partially funded this study. Within the last 5 years, he has received personal fees from GSK, Novartis, Sanofi and Boehringer Ingelheim that are unrelated to the submitted work. To the best of our knowledge, the other authors have no conflict of interest, financial or otherwise.

Author contribution

Liisa Maarit Murray: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Resources; Software; Validation; Visualization; Writing‐original draft; Writing‐review & editing. Gayathri Thillaiyampalam: Data curation; Formal analysis; Investigation; Methodology; Software; Validation; Visualization; Writing‐review & editing. Yang Xi: Conceptualization; Data curation; Methodology; Supervision; Validation. Alexandre S Cristino: Conceptualization; Methodology; Supervision; Writing‐review & editing. John Upham: Conceptualization; Funding acquisition; Investigation; Project administration; Resources; Supervision; Validation; Writing‐review & editing.

Supplementary figures 1‐3

Supplementary tables 1‐7

Click here for additional data file.