Discovering in vivo eQTL interactions with interferon status and drug exposure from a lupus clinical trial

If an expression quantitative trait locus (eQTL) effect is modulated by an environmental stimulus, such as drug exposure or disease status, it can point to key regulatory mediators. In a clinical trial for anti-IL-6 in 157 patients with systemic lupus erythematosus we measured cell counts, interferon (IFN) status, drug exposure and genome-wide gene expression at three time points. First, we confirmed an increase in power using repeat transcriptomic measurements. Then, after detecting 4,976 cis eQTLs, we discovered that 154, 185 and 126 had evidence of significant eQTL interactions with T cell proportion, IFN status and anti-IL-6 drug exposure respectively. Next, we found an enrichment of transcription factor binding motifs interrupted by eQTL interaction SNPs, pointing to regulatory mediators of these environmental stimuli and therefore potential therapeutic targets for autoimmune diseases. In particular, IFN interactions are enriched for IRF1 binding site motifs, while anti-IL-6 interactions are enriched for IRF4 motifs. Finally, we used the drug-eQTL interactions to define an informative drug exposure score, reflecting a drug's effect in an individual patient, thus highlighting the potential for utilizing drug-eQTL interactions in a pharmacogenetic framework.

Here, we predicted that if RNA is queried at multiple time points under different 89 exposure states, repeat measurements could increase power to not only detect eQTLs, 90 but also their interactions with environment. Clinical trials, with their structured study design, may be the ideal setting to detect eQTL 93 interactions with drugs or other environmental perturbations. In clinical trials, it is 94 becoming increasingly common to collect transcriptional and genetic data alongside 95 clinical and physiological data. This extensive phenotyping of the same individuals at 96 multiple time points across conditions can be leveraged to identify eQTL interactions. 97 The profiling of an individual both before and after exposure to a drug provides a unique 98 opportunity to identify in vivo drug-eQTL interactions. In this study, we leverage the power of repeat transcriptional and environmental 117 measurements from a lupus clinical trial to identify eQTL interactions. We discover 118 eQTL interactions with cell count, IFN status and drug exposure. We find the eQTL 119 interaction SNPs are enriched for transcription factor binding motifs (such as IRF1 and 120 IRF4 for IFN and drug exposure respectively) highlighting regulatory mediators of these 121 interactions and potential therapeutic targets.

124
We conducted whole blood high-depth RNA-seq profiling at 0, 12, and 24 weeks in anti-125 IL-6 exposed and unexposed individuals with the Illumina TruSeq protocol. We used a linear mixed model, including repeat measurements with up to three RNA-138 seq assays per patient (Fig. 1A, 379 154 We first tested a type of eQTL interaction that has been examined previously: cell 155 counts 24,25 . We obtained FACS data for 320 samples for which we had 156 contemporaneous RNA-seq profiles (n=152 subjects). We determined the percentage of 157 total lymphocytes that were T cells by gating ( Fig. 2A). We found 154 T cell-eQTL were not inflated, we conducted 1,000 stringent permutations, where we reassigned T 161 cell percentages across samples and retested. This permutation preserved the main 162 eQTL effect, while disrupting interactions that might be present in the data. In no 163 instance did we observe 154 or more interactions at p interact <0.01 (maximum=133), 164 suggesting that the number of observed interactions is highly unlikely to have happened 165 by chance (Fig. 2B, p permute~0 /1,000 = <0.001).  (Fig. 2C, Supplementary Fig. 1). The NOD2 rs1981760 eQTL is an example 170 of an interaction that is dampened by increased T cell count (Fig. 2D, p interact =6.5x10 -5 ), 171 and has separately been shown to vary across cell types 24 183 2). We note that interactions with a proxy gene for IFN status have been described 1 and 184 we find overlap of genes with those reported interactions (Supplementary Fig. 3). To define transcription factors that drive the response to IFN alpha, we sought to 189 identify motifs that explain the differences between magnifier (n=75) and dampener 190 (n=110) eQTLs (Supplementary Fig. 4). We applied HOMER 29 to assess overlap 191 between transcription factor binding motifs and the eQTL interaction SNPs (and SNPs 192 in high linkage disequilibrium (LD, r 2 >0.8) in the cis window) (Methods). We found tight LD (r 2 =0.85 in Europeans) with rs6494127, which interrupts the GAAA core of the 199 IRF1 motif (Fig. 3C), and likely disrupts IRF1 binding 30 . We observe greater expression 200 of GTF2A2 in individuals with the rs2306355 A allele compared to G; this difference is 201 dampened in IFN low individuals (Fig. 3D).

203
Discovery of eQTL interactions with drug exposure 204 We then examined whether IL-6 blockade alters the relationship between genomic 205 variation and gene expression and induces drug-eQTL interactions. We observed 126 9 drug-eQTL interactions with p interact <0.01 (Supplementary Table 1). Following the same 207 permutation strategy as above, we found a median of 77 interactions with p interact <0.01 208 (maximum=117) from 1,000 permutations. This suggests that about half of our drug-209 eQTL interactions likely represent real biological phenomena, and not statistical artifact 210 (Supplementary Fig. 5). These drug-eQTL interactions showed little overlap with the 211 interactions observed for T cell count or IFN status (Supplementary Fig. 6). We note 212 biologically relevant drug-eQTL interactions for IL10 (Supplementary Fig. 7), an anti-  Fig. 4C, p=3.2x10 -7 , binomial test). We were concerned that free 250 IL-6 and drug exposure are not independent in this dataset, and that this concordance 251 might be in part due to the connection between free IL-6 protein levels and IL-6 252 blockade (Supplementary Fig. 10). To assess whether free IL-6 offers independent interaction effects that were consistent after accounting for drug effect, we first modeled 254 free IL-6 levels to account for the presence or absence of drug, and then we assessed 255 interactions with residual IL-6 levels. Again, we observed a significant number of 256 interactions in a consistent direction with residual IL-6 levels, which are independent of 257 drug exposure (p=0.03, Supplementary Fig. 11).

259
Drug score for assessing drug exposure 260 We speculated that these drug-eQTL interactions could be used in a clinical 261 pharmacogenetic context to assess effective drug exposure for patients. We defined a 262 simple drug exposure score using the 126 drug-eQTL interactions (Methods). For each 263 RNA-seq sample, we assessed whether the expression of the interaction gene was 264 more consistent with the drug exposed or unexposed state for the corresponding 265 interaction SNP genotype. Samples more consistent with the drug-exposed state are 266 assigned a larger drug exposure score. Unsurprisingly, we found a difference in drug 267 exposure score between the unexposed and exposed samples (Supplementary Fig.   268 12) (r s =0.79, p=6.9x10 -81 ); these differences reflect the fact that the eQTLs were 269 themselves identified by examining samples with and without drug exposure. However, 270 while we did not utilize the administered drug dose to identify drug-eQTL interactions, 271 we found a significant correlation between drug dose (10, 50 or 200mg) and drug 272 exposure score (r s =0.16, p=0.018) in the drug-exposed samples (Fig. 4D).  Fig. 7). Previous studies have highlighted a role 314 for IRF4 in the pathogenesis of autoimmune diseases in mouse and humans. For 315 example in a murine model of SLE, IRF4 knockout mice did not develop lupus 316 nephritis 33 . In humans, IRF4 is associated with RA 34 , a disease in which anti-IL-6 317 treatment has been successful 21 . Our findings provide further support that IRF4 could 318 be a potential therapeutic target for autoimmune diseases such as RA where anti-IL-6 is 319 effective 35 . As the regulatory genome continues to be mapped 36,37 and available binding 320 sites for different regulator proteins are specifically defined, these approaches to define 321 environmental mechanisms will become more potent. For example, it will become easier to connect groups of interaction eQTLs being driven by specific transcription factors or a 323 subset of regulatory elements. Finally, we argue that drug-eQTL interactions can augment pharmacogenetic strategies 339 and may be informative for patient response. For many biologic medications, predictive 340 pharmacogenetics has been challenging; for example, studies to define genetic or non-341 genetic biomarkers of anti-TNF response have not been successful 39,40 . Our drug 342 exposure score, based on a novel approach using SNP-gene pairs from drug-eQTL 343 interactions, might reflect the biological activity that a medication is having upon an 344 individual, and may be modeling an effective medication activity level. This score may therefore be helpful in stratifying individuals when assessing response to a medication 346 for example those with a higher drug exposure score may have a better response to 347 treatment. We note a limitation of this study is that the drug itself did not achieve its 348 primary efficacy endpoint of improving SLE outcomes. Hence, while the drug exposure 349 score for this study tracked with the biological effect of the drug (reducing free IL-6 350 protein levels), it might not be useful for SLE specifically. However, such a scoring   Hardy-Weinberg Equilibrium (p < 1x10 -7 ), had a minor allele frequency <5%, 393 missingness >2% or a heterozygosity rate greater than 3 standard deviations from the 394 mean (PLINK 44,45 ). For mapping eQTLs, we removed SNPs on the Y chromosome. 395 Following QC, we used 608,017 variants for further analysis. We removed one sample 396 with high missingness and outlying heterozygosity rate from further analysis.

398
Cell counts 399 We collected blood samples for cytometry analysis at weeks 0, 12 and 24. Samples

Drug exposure
Samples were assigned as unexposed (placebo or week 0 samples) or drug exposed 413 (week 12 and week 24 samples in the drug groups).

415
Free IL-6 protein levels 416 We determined free IL-6 protein levels from serum using a commercial sandwich ELISA 417 selected for binding only free IL-6. The assay was validated according to FDA   474 In the SLE cohort, we classified 4,976 cis eQTL genes (p<2.3x10 -8 ). The z-score for the 475 most associated SNP for each of these genes was compared to the z-score from a 476 previously published eQTL dataset from whole blood from 2,166 healthy individuals 1 . 477 4,250/4976 SNP-gene pairs (85.4%) were also reported in the BIOS dataset Conditional analysis for IL-6 protein levels 494 We modeled the relationship between free IL-6 protein levels and drug exposure using 495 a linear model. We used the residuals from this model in our interaction linear mixed 496 model to identify IL-6 protein interactions independent of drug exposure.

498
Drug exposure score 499 We used linear discriminant analysis to assign a drug exposure score for each sample. 500 A score was calculated for each gene (see equation below) and then the final drug 501 exposure score is the average across the 126 drug-eQTL genes.
Where G is gene expression for a given sample, G Unexp is predicted mean gene