Host Genetic Variation Influences Gene Expression Response to Rhinovirus Infection

Summary

Rhinovirus (RV) is the predominant cause of the common cold. However, infections with RV result in a broad spectrum of effects ranging from asymptomatic infections to severe lower respiratory illnesses. We hypothesized that diversity in response to RV-infections is, at least in part, due to variation in the host genome. To address this, we mapped the genetic variations that are associated with gene expression response (reQTLs) to RV-infection in PBMCs. Here, we report local reQTLs for 38 genes including those with known functions in viral response such as UBA7, OAS1, IRF5 and those that have been previously associated with immune and RV-related diseases (e.g., ITGA2, MSR1, GSTM3). We also show that reQTL regions are enriched for binding sites of the virus-activated STAT2 transcription factor, suggesting a potential mechanism of action for five of the reQTLs identified. Overall, the reQTLs we identified represent promising candidates to affect individual’s immune response to RV infections and further targeted studies of the reQTL regions might lead to improved control and treatment of RV-associated immune and 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).

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).

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 (r²>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).

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, 4x10⁶ PBMCs were treated with media alone for 24 hours and 4x10⁶ 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 log₂-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; r²>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 r² 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.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18