Host genetic variation influences gene expression response to rhinovirus infection.

Rhinovirus (RV) is the most prevalent human respiratory virus and is responsible for at least half of all common colds. RV infections may result in a broad spectrum of effects that range from asymptomatic infections to severe lower respiratory illnesses. The basis for inter-individual variation in the response to RV infection is not well understood. In this study, we explored whether host genetic variation is associated with variation in gene expression response to RV infections between individuals. To do so, we obtained genome-wide genotype and gene expression data in uninfected and RV-infected peripheral blood mononuclear cells (PBMCs) from 98 individuals. We mapped local and distant genetic variation that is associated with inter-individual differences in gene expression levels (eQTLs) in both uninfected and RV-infected cells. We focused specifically on response eQTLs (reQTLs), namely, genetic associations with inter-individual variation in gene expression response to RV infection. We identified local reQTLs for 38 genes, including genes with known functions in viral response (UBA7, OAS1, IRF5) and genes that have been associated with immune and RV-related diseases (e.g., ITGA2, MSR1, GSTM3). The putative regulatory regions of genes with reQTLs were enriched for binding sites of virus-activated STAT2, highlighting the role of condition-specific transcription factors in genotype-by-environment interactions. Overall, we suggest that the 38 loci associated with inter-individual variation in gene expression response to RV-infection represent promising candidates for affecting immune and RV-related respiratory diseases.

Introduction

Rhinovirus (RV) is the most prevalent human respiratory pathogen [1]. It was discovered as the predominant cause of the common cold over 50 years ago [2]. Longitudinal studies indicate that nearly all individuals experience at least one RV infection by two years of age [3]. Each year following, pre-school age children experience six [4] and adults experience two to three [5] RV infections on average. Recent studies have shown that infection with RV results in a broad spectrum of illness severity, ranging from asymptomatic infections to severe lower respiratory illnesses such as bronchiolitis and pneumonia [6]. In addition, RV infections contribute to the morbidity of chronic respiratory illnesses such as asthma, chronic obstructive pulmonary disease (COPD), and cystic fibrosis (CF) [6]. The diversity in response to RV infection is likely attributable to, at least in part, inter-individual variation in the host genome. Indeed, genetic variation in the promoter region of the IL-10 gene was shown to influence severity of RV illnesses in a small sample of 18 subjects [7] and polymorphisms at the 17q12-q21 asthma locus were associated with both the occurrence and number of RV wheezing illnesses in early life [8]. Beyond these few associations however, the genetic and/or mechanistic basis for the vast inter-individual variation in the response to RV infection is not well understood.

Many cell types are involved in the immune response to RV [9]. Genome-wide gene expression response to RV infection was previously studied in bronchial epithelial cells [10,11] and in nasal epithelial scrapings [12]. While the nasal mucosa is considered the primary site of RV replication [13], RV genomes were found in pericardial fluid, stool, urine, plasma, and serum samples of children with respiratory illnesses, suggesting that systemic infection of RV occurs [14-16]. In addition, RV infection induces cytokine production from monocytes and macrophages without productive viral replication [17], raising the possibility that RV-associated respiratory illnesses may result from virus-induced inflammatory cytokines rather than to cytopathic effects of RV per se [18]. To date, however, there have been no genome-wide studies of gene expression response to RV infection in peripheral blood mononuclear cells (PBMCs), or studies characterizing the genetic architecture of inter-individual regulatory variation in gene expression response to RV.

To begin addressing this gap, we have collected and analyzed gene expression data in PBMCs of 98 individuals, before and after RV infection in vitro. Our study design allowed us to provide a comprehensive view of the regulatory variation involved in gene expression levels (eQTLs) in uninfected and RV-infected PBMCs. We also identified genetic variations that are specifically associated with differences in gene expression response (reQTLs) to RV infection; these are loci that interact, directly or indirectly, with the infection process and likely include a subset of loci that contribute to the inter-individual variation in the clinical response to RV infection.

Results

Rhinovirus infection has a systematic effect on gene expression

We obtained genome-wide array-based gene expression data from uninfected and RV-infected PBMCs from 98 unrelated adults (GEO accession number: GSE53543). We detected as expressed 13,881 autosomal probes (targeting 10,893 genes; see Materials and Methods). Across the 196 samples (uninfected and RV-infected paired samples from 98 individuals), gene expression profiles clearly clustered into two major groups by treatment status (by PCA; see Methods; S1 Fig and S2 Fig). We next identified the differentially expressed genes between uninfected and RV-infected PBMCs; of the 10,893 expressed genes, 2,242 were up-regulated and 3,779 were down-regulated in RV-infected compared to uninfected PBMCs (at the Bonferroni corrected significance threshold of P<4.60x10-6; Fig. 1A, S1 Dataset). We note that with a samples size of 98 individuals we have considerable power to detect differences in gene expression levels. Indeed, the vast majority of effect sizes we detected as significant were minor (Fig. 1A).

Fig 1

Identification and functional characterization of RV-responsive genes in PBMCs.

(A) Volcano plot showing the Log2 Fold change (x-axis) and the—Log10 P value (y-axis) of gene expression between uninfected and RV-infected PBMCs. Genes that were not classified as differentially expressed are shown in black. Genes that were significantly differentially expressed (P<4.6x10-6) but displayed <2-fold change are shown in blue; genes that were both significantly differentially expressed and had ≥2-fold change are shown in red. (B) Gene ontology (GO) enrichment results of the genes that were both significantly differentially expressed and had ≥2-fold change (see S3 Table for results including all GO terms). (C) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment results of the genes that were both significantly differentially expressed and had ≥2-fold change (see S4 Table for results including all KEGG terms).

We focused, therefore, on the 271 genes that were both significantly differentially expressed and showed at least two-fold difference in expression between uninfected and RV-infected PBMCs. Although a two-fold threshold is arbitrary, it is consistent with the thresholds used in other studies of gene expression response to RV [11,12] and thus provides some measure of consistency with previous reports. In fact, considering our data in the context of other studies we observed a large overlap between the genes that respond to RV in different cell types (86.5% overlap between PBMCs and bronchial epithelial cells [S1 Table] and 73.3% overlap between PBMCs and nasal epithelial scrapings [S2 Table]). Additionally, pathway analysis of the 271 genes revealed enrichment for immune and viral response related pathways, as expected (Fig. 1B, Fig. 1C, S3 Table, S4 Table).

Mapping eQTLs in uninfected and rhinovirus-infected PBMCs

We next identified local (putative cis) eQTLs in uninfected and RV-infected PBMCs by testing for associations between expression levels and genetic variation at loci within 1 Mb windows from the nearest annotated end of each expressed gene. We performed this analysis separately in uninfected and RV-infected PBMCs, and in both cases we used the SNP with the minimum P value observed for each gene to assess the evidence of an eQTL at that locus. At a false discovery rate (FDR) of 5% based on permutations (which also consider the minimum P value for each gene), we identified 521 genes with local eQTLs in uninfected PBMCs (Fig. 2A, S2 Dataset) and 523 genes with local eQTLs in RV-infected PBMCs (Fig. 2A, S2 Dataset).

Fig 2

Identification of local eQTLs and reQTLs.

(A) Manhattan plots of local eQTLs in uninfected cells (in blue), local eQTLs in RV-infected cells (in red), and local reQTLs (in green). Most significant local eQTL or local reQTL P value for each gene (y-axis) is displayed in the order of chromosomal positions of the genes detected as expressed in our study (x-axis). (B) Examples of genes with reQTLs. In each plot, genotype at the reQTL is shown on the x-axis and expression level (Log2 Expression) or response in gene expression (Log2 Fold Change) is shown on the y-axis. Sample sizes for each genotype group are shown in parentheses under the x-axis.

For each significant local eQTL-gene expression pair in uninfected and in RV-infected PBMCs, we compared the evidence of eQTL association across conditions (S3 Fig). Pearson correlation of eQTL association P values was 0.84 for the 521 eQTLs in uninfected cells (P<10–15) and 0.81 for the 523 eQTLs in RV-infected cells (P<10–15), suggesting that a significant proportion of genetic regulation of gene expression is maintained between uninfected and RV-infected PBMCs.

Mapping response eQTLs (reQTLs)

We reasoned that SNPs that are specifically associated with gene expression response to RV infection (reQTLs) are likely to have a role in RV-specific response. To identify reQTLs, we tested associations between the expression response (estimated as Log2 fold change in gene expression response to RV infection) and genetic variation in loci that are within 1 Mb of the nearest end of each expressed gene. This analysis revealed local reQTLs for 38 genes at the FDR 5% threshold (Fig. 2A, S2 Dataset). For 25 of the genes with reQTLs, genotypic effect on gene expression at FDR 5% threshold was present only in uninfected or only in RV-infected cells (S2 Dataset). For the remaining 13 genes, genotypic effect on gene expression was present in both uninfected and RV-infected cells at FDR 5% threshold, but the effect size of association was significantly different between conditions (S2 Dataset). Overall, expression levels of genes with reQTLs were slightly higher in the condition where absolute genotypic effect size of the reQTL was larger (S4 Fig). This analysis raises the possibility that slight differences in power to map eQTLs may explain a subset of our observations, but we note that the average difference in expression level is small and is not consistent across all genes with reQTLs.

The 38 genes with reQTLs included those with known functions in viral response, such as UBA7, OAS1, IRF5 (Fig. 2B). The protein product of UBA7 (ubiquitin-like modifier activating enzyme 7) is involved in the activation of a critical antiviral protein ISG15 (interferon-stimulated gene 15) [19,20]. OAS1 (2’-5’-oligoadenylate synthase 1) activates latent RNase L following viral infections and results in degradation of viral RNA [21,22]. Similarly, IRF5 (interferon regulatory factor 5) protein is a direct transducer of virus-mediated signaling [23]. In addition, a subset of the 38 genes with local reQTLs have previously been associated with immune or RV-associated diseases (such as asthma, COPD, and CF; as catalogued by the Genetic Association Database [24]). Specifically, eight genes (IRF5, UBA7, TMTC1, GSTM3, FBN2, OAS1, PODXL, ITGA2) were previously associated with immune system diseases in at least one study [25-29]. Notably, three genes (IRF5, GSTM3, FBN2) have been associated with asthma [29-31], two (GSTM3, MSR1) with COPD [32,33] and one (GSTM3) with CF [34]. It should be noted, however, that while these observations represent a slight enrichment compared to genome-wide expectations (S5 Fig), enrichment P values were not significant.

In an attempt to fine map the causal reQTLs, we next repeated reQTL mapping for the 38 genes with significant reQTLs using imputed genotype data (see Material and Methods for details of imputation). We compared the most significant reQTL for each gene based on imputed vs. genotyped data (S5 Table). For 16 of the genes, the most significant reQTL based on genotyped data remained the most significant reQTL based on imputed genotype data. However, for 12 of the genes, multiple SNPs had the smallest reQTL P value in the imputed genotypes, suggesting that fine mapping efforts using imputed genotype data in a sample size of 98 individuals is not sufficient to identify the causal reQTL.

Enrichment of STAT2 binding sites among reQTL loci

Previous studies have suggested that reQTLs can often be found within the binding sites of transcription factors that have different levels of activity between the conditions being tested [35,36]. We examined this in our data by considering the overlap between known protein binding sites (based on ENCODE ChIP-Seq data [37]) and the 38 reQTLs identified in our study. This analysis is challenging because of the uncertainty in identifying the causal reQTL association. We therefore considered for this analysis the reQTLs and all other common SNPs in high LD with the reQTLs (r2>0.8 and MAF>0.1 based on 1000 Genomes Phase I European Population) as input SNPs. We then extracted 1,000 control SNPs (matched based on MAF, gene density, distance to nearest gene, number of SNPs in LD) per each input SNP and examined the overlap between protein binding sites and each unique SNP in the input and control SNP lists. This analysis revealed that STAT2 was the most substantially enriched transcription factor (16.30 fold; P<10–6; S6 Table and Fig. 3A) with binding sites overlapping reQTL loci SNPs relative to the control SNPs (see Materials and Methods for details of enrichment analysis).

Fig 3

STAT2 binding sites are enriched among reQTL loci.

(A) Observed and expected numbers of variants that lie within STAT2 binding sites i) across all human cell types/conditions available in ENCODE ChIP-Seq data, ii) in IFNα-30 minute treated human K562 cell line, iii) in IFNα-6 hour treated human K562 cell line. (B) Regional eQTL association plots of SLFN5 in uninfected and RV-infected cells. See S6 Fig for regional eQTL association plots of the remaining four genes with significant reQTL loci SNPs that reside in STAT2 binding sites. In all five cases, the eQTL association was significant only in RV-infected cells.

STAT2 is a critical regulator of anti-viral response [38]. In our study, STAT2 gene expression increased by 3.68 fold (S1 Dataset; P<10–77) in response to RV infection. Based on this observation, we hypothesized that reQTL loci SNPs that lie within STAT2 binding sites should have stronger regulatory effects on gene expression in RV-infected cells relative to uninfected cells. As expected, for all five targets of STAT2 among our reQTL regions, the eQTL effects on gene expression were stronger in RV-infected cells compared to uninfected cells (Fig. 3B and S6 Fig). Consistent with this observation is that all the STAT2 ChIP-Seq signals (based on ENCODE ChIP-Seq data) in our reQTL loci were identified in the IFNα treated human K562 cell line (S7 Table and Fig. 3A). Additionally, our results imply that the SNPs residing in STAT2 binding sites are more likely to be the causal SNP in the reQTL loci of EXOSC9, SLFN5, PRR24, OAS1, and ARL5B. In fact, causality of the reQTL of SLFN5 gene, rs11080327, was recently demonstrated by luciferase reporter assays [35].

Little evidence for trans reQTLs

Lastly, we searched for distant (putatively trans) eQTLs and reQTLs by testing associations between SNP-gene combinations for which the SNP distance from the nearest end of the gene was more than 5 Mb. We used the minimum P value observed per gene to assess the significance of the distant eQTL or reQTL. We identified 11 and 15 genes with distant eQTLs in uninfected and in RV-infected PBMCs, respectively (at an FDR 5% based on permutations; S3 Dataset). We found no direct evidence for distant reQTLs when we considered genetic associations with the expression response to infection (S3 Dataset), yet 6 of the trans eQTLs were specific to either the infected or uninfected cells. Because trans eQTL findings based on microarray expression data tend to suffer more from a high degree of false-positives, partly due to cross-hybridizations, we considered carefully each significant finding. For each gene that was putatively associated with a trans eQTL, we re-mapped each of its probes within 2 Mb of the trans eQTL using SHRiMP [39], as previously described [40]. We discarded all trans eQTLs whose associated trans-probe mapped in its vicinity. After exclusions, only 4 and 6 genes with distant eQTLs in uninfected and RV-infected cells, respectively, remained; 3 of them were common to both conditions (S8 Table). Thus, 4 trans eQTLs may be reQTLs, yet considerations of incomplete power for detecting trans eQTLs in our study call for caution in this interpretation.

Discussion

Mapping of condition-specific or response eQTLs in a genome-wide level is an emerging area of research with the promise of understanding biology of infectious diseases [41,42], response to pharmaceutical treatment [36], and inter-individual variation in immune response more broadly [35,43]. Here, we report 38 local reQTLs that are associated with inter-individual variation in gene expression response to RV infection. These loci are likely to interact, directly or indirectly, with the infection process. Because infection with rhinovirus causes significant changes in regulatory activity of the identified reQTLs, these variants represent promising candidates for susceptibility and response to RV infections in vivo. In support of this reasoning, we pointed to eight genes with RV-reQTLs that have been previously associated with immune diseases, and four that were previously associated with respiratory diseases (either asthma, COPD, or CF).

We note that while the numbers of disease-associated genes among genes with reQTLs were greater than the genome-wide expectations, none of the enrichment P values were significant. That said, it is important to emphasize that RV infections, in general, are thought to affect the morbidity of the chronic respiratory illnesses rather than the risk of developing the disease [6]. Therefore, it is possible that the reQTLs identified here are more likely to influence the severity of the respiratory illnesses rather than the occurrence, which has been the focus of the majority of the association studies performed thus far. Similarly, risk of developing the common cold might be more directly influenced by the reQTLs identified here because RV is causally associated with development of the common cold. However, to our knowledge, there have been no association studies on the frequency or severity of common colds.

It is also possible that at least some of the reQTLs identified here are functional in response to a broader range of stimuli, including other pathogens. In fact, six (SLFN5, ARL5B, SPTLC2, IRF5, ADCY3, CCDC146) of the 38 genes implicated in our reQTL mapping were also identified in a recent reQTL mapping study for Escherichia coli lipopolysaccharide, influenza, and interferon-β in dendritic cells [35]. Further comparative studies of reQTLs will be necessary to disentangle stimulus-specific and shared reQTLs. We also note that some of the reQTLs identified in our study may be confounded due to cell type heterogeneity in PBMCs. In our study, we were unable to count cell subsets of PBMCs before and after RV infection. If for instance, cell subset proportions of PBMCs change in response to RV infection, statistical power to identify cell-type specific eQTLs may differ between uninfected and RV-infected PBMCs and this may potentially lead to identification of cell-type specific eQTLs as reQTLs. Similarly, it is possible that cell type heterogeneity in PBMCs may mask identification of cell-type specific reQTLs, especially those that are specific to rare cell subsets of PBMCs. However, we also appreciate the fact that gene expression response in a complex system of interacting cells as in PBMCs might be more relevant to true physiological responses than those observed in purified cell subsets.

The enrichment of STAT2 binding sites among reQTL regions highlights the role of condition-specific transcription factors in gene-by-environment interactions. Our results suggest that transcription factor activation upon RV-infection reveals SNPs with regulatory activity that could not be identified in uninfected PBMCs. This hypothesis is also supported by the fact that all STAT2 ChIP-Seq signals in our reQTL regions were identified in the IFNα treated human K562 cell line. IFNα is involved in the innate immune response against viral infections and our results therefore suggest that RV infection may activate STAT2 through the IFNα signaling pathway.

In conclusion, we have provided a comprehensive genome-wide view of host genetic variation that is associated with gene expression response to rhinovirus infection. The reQTLs identified here are promising candidates to influence both the frequency and the severity of RV related respiratory illnesses. Additionally, our results contribute to the field of genotype-by-environment interactions and might further help to disentangle stimulus-specific and shared reQTLs.

Materials and Methods

Ethics statement

One hundred unrelated adult volunteers (49 males and 51 females; age range 19–60) were recruited between July and November 2011 to study the genotype-specific effects of RV infection on gene expression patterns in PBMCs. Informed written consent was obtained from each study participant. This study was approved by the Institutional Review Board at the University of Chicago.

Sample collection and experimental design

Twenty ml of blood was drawn from each participant. PBMCs were isolated from whole blood samples by Ficoll-Paque separation [44]. From each subject, 4x106 PBMCs were treated with media alone for 24 hours and 4x106 PBMCs were treated with media containing RV16 for 24 hours. The multiplicity of infection was 10 plaque-forming units per cell.

DNA extraction and genome-wide genotyping

DNA was extracted on the day of sample collection using QIAamp DNA Blood Mini Kit; concentrations were measured on a Nanodrop ND-100 Spectrophotometer. Genotyping of 100 individuals was performed using the Axiom Genome Wide Human Array Plate–CEU, which interrogates 669,059 SNPs. Genotype calls were extracted from the raw data using Affymetrix Power Tools software. One individual with genotype call rate less than 97% and one individual that failed Affymetrix gender call were excluded from all further analyses. After exclusions, the sample size was 98 (49 males and 49 females; age range 19–60). Following quality-control checks (Hardy-Weinberg equilibrium P>10–6, MAF>0.1, SNP call rate>95%), 382,855 SNPs were retained and 373,312 autosomal SNPs with unique SNP identifiers were used in subsequent analyses.

Global proportions of European, Asian, and African ancestry in our samples were estimated by using the program ADMIXTURE [45] and assuming 3 ancestral populations (K = 3) (S7 Fig). Subjects from the Phase 3 HapMap CEU (Utah residents with Northern and Western European ancestry from the CEPH collection), CHB (Han Chinese in Beijing, China), JPT (Japanese in Tokyo, Japan), and YRI (Yoruba in Ibadan, Nigeria) were included as reference populations. 85 of the individuals had over 80% European ancestry, six of them had over 80% Asian ancestry and the remaining seven individuals had relatively more mixed ancestry fractions. Ancestry estimates were taken into account in further analyses to correct for the potential effects of ancestry on gene expression profiles.

Genotype imputation

373,312 SNPs that passed the genotyping quality-control checks were used to perform the pre-phasing of the chromosomes using SHAPEIT (v2.r790) [46]. Imputation was performed using IMPUTE2 (2.3.1) [47] over genomic regions of 5 Mb, as recommended. For both pre-phasing and imputation, 1000 Genomes Phase 3 data were used as the reference panel.

78,091,231 autosomal variant sites were imputed across 98 individuals. After quality-control checks (Hardy-Weinberg equilibrium P>10–6, MAF>0.1, SNP call rate>95%), 3,722,989 of the imputed variants were retained.

RNA extraction and gene expression profiling

Total RNA was extracted after 24-hour incubation, using the RNeasy Plus Mini Kit; concentrations were measured on a Nanodrop ND-100 Spectrophotometer and quality was assessed using an Agilent 2100 Bioanalyzer. Genome wide gene expression profiling of uninfected and RV-infected PBMCs was obtained using Illumina HumanHT-12 v4 Expression BeadChip arrays, which targets 47,305 probes. The cDNA synthesis, labeling, and hybridization of RNA to the microarrays were performed at the University of Chicago Functional Genomics Core.

The process of quality-control checks (probes that mapped to unique Ensembl gene IDs, probes that did not contain any HapMap SNPs with MAF>0.01 in the CEU population, probes that targeted autosomal chromosomes) resulted in retention of 26,440 probes. Among those, 13,881 probes that were detected as expressed in PBMCs (detection P<0.05 in at least 25% of the samples) were used in subsequent analyses. Low-level microarray analyses were performed in R (http://www.R-project.org), using the Bioconductor software package lumi [48]. Probe intensity estimates were log2-transformed and rank-invariate normalized.

Linear models were used to test the relationship between each known covariate (gender, virus batch, processing day, chip number, PBMC count, age, and ADMIXTURE’s [45] Q estimates (the ancestry fractions)) and the principal components that explain at least 5% of the total variance in the gene expression data. (S1 Fig). Processing day was significantly associated with the Principal Component 2 and hence it was included as an “adjustment variable” when performing SVA [49]. SVA analysis yielded no significant surrogate variables when “adjustment variable” of processing day was used with “variable of interest” of treatment. The effects of processing day were regressed out of the gene expression data prior to further analyses (S2 Fig).

Median probe intensity estimates per gene were calculated as the expression levels for 10,893 genes and used in all further analyses.

Identification of differentially expressed genes between uninfected and RV-infected cells

Differentially expressed genes between uninfected PBMCs and RV-infected PBMCs were identified using a paired t-test in R statistical environment. Significance was calculated using the Bonferroni correction at α = 0.05 (P<4.60x10-6).

Gene ontology, KEGG pathway enrichment analyses

The DAVID bioinformatics database [50] was used to test for enrichment of Gene Ontology (GO) categories (BP_ALL) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways among the genes that were both statistically differentially expressed and had ≥2-fold change in response to RV infection. 10,893 genes that were detected as expressed in our data were used as the background set in all enrichment analyses. Enrichment P values were calculated using a modified Fisher Exact test (EASE Score).

Genetic mapping of gene expression and gene expression response

Associations between local (defined as SNP-gene pairs within 1 Mb window from the nearest end of the gene) and distant (defined as SNP-gene pairs for which the SNP distance from the nearest end of the gene was more than 5 Mb) SNPs and gene expression in uninfected and in RV-infected cells were tested using linear regressions with additive genotype effects and taking ancestry estimates into account. Similarly, associations between local and distant SNPs and gene expression response (Log2 fold change in gene expression response to RV infection) were tested using linear regressions with additive genotype effects and taking ancestry estimates into account. All the analyses were performed as implemented in “linear” model of the Matrix eQTL package [51]. The minimum P value observed for each gene was recorded and used as the evidence of eQTL or reQTL association.

To estimate the FDR, each of the phenotype data (gene expression or gene expression response) were shuffled 100 times and linear regressions were repeated using each set of permuted data. The minimum P value for each gene was recorded and used as the empirical null distribution. Permutation-based FDR was calculated using fdrci package [52] in R.

Analysis of gene-disease associations

Genes previously reported to be associated with immune diseases (Disease Class) and RV-related respiratory illnesses (asthma, COPD, CF) (Disease Terms) were downloaded from the Genetic Association Database (GADCDC data as of 08/18/2014). Enrichment P values were calculated using Fisher's exact test.

ChIP-Seq enrichment analysis

SNP ‒ human ChIP-Seq (based on ENCODE data [37]) annotation data were downloaded from HaploReg v2 [53] in December 2014. For each ChIP annotation across all available cell types and conditions, the frequency of overlap between unique reQTL loci SNPs was compared with the frequency of overlap with a set of unique control SNPs.

reQTL loci SNPs included the list of reQTLs and all other common SNPs in high LD with the reQTLs; r2>0.8 and MAF>0.1 based on 1000 Genomes Phase I European Population. 1,000 control SNPs (matched based on MAF [±5% point], gene density [±50%], distance to nearest gene [±50%], number of SNPs in LD [±50% using r2 0.5]) per each input SNP were pulled using the SNPsnap webserver [54]. Binomial P values and fold-enrichments were calculated for proteins with at least two ChIP-Seq signal overlapping with the reQTL loci SNPs.

For STAT2 ChIP annotation, the frequency of overlap between reQTL loci SNPs and the matching control SNPs was additionally compared within IFNα-30 minute and IFNα-6 hour treated human K562 cell line.

Quality control check of trans eQTLs

To minimize the false-positive trans eQTL findings, each probe targeting the genes with significant trans eQTLs was re-mapped using SHRIMP [39] with relaxed mapping setting that was previously described; match score of 10, mismatch score of 0, gap open penalty of −250, gap extension penalty of −100, and minimal Smith-Waterman score of 30% [40]. Trans eQTLs whose associated trans-probe mapped in its vicinity were excluded.

Clustering of gene expression data.

(A) Heatmap clustering, and (B) Principal component plot of 13,881 probes after log2-transformation and rank-invariate normalization. (C) P values from linear models testing the relationship between each known variable (potential covariates and variable of interest) and the principal components that explain at least 5% of the total variance in the gene expression data are shown. P values that are significant (after Bonferroni correction at α = 0.05) are highlighted in red.

(PDF)

Click here for additional data file.

Clustering of gene expression data after regressing out processing day.

(A) Heatmap clustering, and (B) Principal component plot of 13,881 probes after log2-transformation, rank-invariate normalization, and regressing out processing day. (C) P values from linear models testing the relationship between each known variable (potential covariates and variable of interest) and the principal components that explain at least 5% of the total variance in the gene expression data are shown. P values that are significant (after Bonferroni correction at α = 0.05) are highlighted in red.

(PDF)

Click here for additional data file.

Comparisons of the evidence of eQTL association across conditions.

(A) For 521 significant local eQTL-gene expression pair in uninfected cells, P values in uninfected cells are shown on the x-axis and P values in RV-infected cells are shown on the y-axis. Pearson correlation of eQTL association P values was 0.84 (P<2.2x10-16). (B) For 523 significant local eQTL-gene expression pair in RV-infected cells, P values in RV-infected cells are shown on the x-axis and P values in uninfected cells are shown on the y-axis. Pearson correlation of eQTL association P values was 0.81 (P<2.2x10-16). (C) For 521 significant local eQTL-gene expression pair in uninfected cells, Beta (effect size estimate) values in uninfected cells are shown on the x-axis and Beta values in RV-infected cells are shown on the y-axis. Pearson correlation of Beta values was 0.94 (P<2.2x10-16). (D) For 523 significant local eQTL-gene expression pair in RV-infected cells, Beta values in RV-infected cells are shown on the x-axis and Beta values in uninfected cells are shown on the y-axis. Pearson correlation of Beta values was 0.93 (P<2.2x10-16).

(PDF)

Click here for additional data file.

Expression levels of genes with reQTLs in the condition where absolute genotypic effect size of the reQTL was larger vs. in the condition where absolute genotypic effect size of the reQTL was smaller.

(PDF)

Click here for additional data file.

Percentage of disease associated genes.

Percentage of genes associated with (A) Immune diseases (B) Asthma (C) Chronic Obstructive Pulmonary Disease (COPD) (D) Cystic Fibrosis (CF). In each panel, ‘all genes’ category (in black) refers to 10,893 genes detected as expressed and ‘genes with reQTLs’ category (in green) refers to 38 genes that had a significant local reQTL in our study. Enrichment P values (Fisher’s exact test) were 0.17, 0.17, 0.23, and 0.12, respectively.

(PDF)

Click here for additional data file.

Regional eQTL association plots of genes with reQTL loci variants that reside in STAT2 binding sites.

(PDF)

Click here for additional data file.

Global proportions of European, Asian, and African ancestry in our study subjects were estimated using ADMIXTURE, assuming 3 ancestral populations.

Subjects from the Phase 3 HapMap CEU (Utah residents with Northern and Western European ancestry from the CEPH collection), CHB (Han Chinese in Beijing, China), JPT (Japanese in Tokyo, Japan), and YRI (Yoruba in Ibadan, Nigeria) were included as reference.

(PDF)

Click here for additional data file.

Overlap between RV-responsive genes in bronchial epithelial cells (BECs) and in PBMCs.

Among the genes that were targeted in both studies, 86.5% (32 out of 37) with ≥2-fold increase in response to RV infection in BECs also showed ≥2-fold increase in response to RV infection in PBMCs.

(PDF)

Click here for additional data file.

Overlap between RV-responsive genes in nasal epithelial scrapings (NESs) and in PBMCs.

Among the genes that were targeted in both studies, 73.3% (22 out of 30) with ≥2-fold increase in response to RV infection in NESs also showed ≥2-fold increase in response to RV infection in PBMCs.

(PDF)

Click here for additional data file.

Gene ontology (GO) enrichment analysis for the 271 genes that were both statistically differentially expressed and had a fold change of ≥2 in response to RV infection.

Enriched GO categories at P value cut-off of 10–5 are shown.

(PDF)

Click here for additional data file.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the 271 genes that were both statistically differentially expressed and had a fold change of ≥2 in response to RV infection.

Enriched KEGG pathways at P value cut-off of 0.05 are shown.

(PDF)

Click here for additional data file.

Comparison of most significant reQTL based on imputed and genotyped data.

(PDF)

Click here for additional data file.

ChIP-Seq enrichment results of input SNPs relative to the background SNPs based on all available human cell types and treatment conditions in the ENCODE data.

(PDF)

Click here for additional data file.

reQTL loci variants that overlap with STAT2 ChIP-Seq annotation.

LD (r2) is based on 1000 Genomes Phase 1 EUR population.

(PDF)

Click here for additional data file.

Quality control check of genome-wide significant distant eQTLs (A) in uninfected cells (B) in RV-infected cells.

(PDF)

Click here for additional data file.

Summary of the results of differential expression analysis.

(XLSX)

Click here for additional data file.

Summary of the results of local eQTL and reQTL mapping.

(XLSX)

Click here for additional data file.

Summary of the results of distant eQTL and reQTL mapping.

(XLSX)

Click here for additional data file.