Whole Genome DNA and RNA Sequencing of Whole Blood Elucidates the Genetic Architecture of Gene Expression Underlying a Wide Range of Diseases

To create a scientific resource of expression quantitative trail loci (eQTL), we conducted a genome-wide association study (GWAS) using genotypes obtained from whole genome sequencing (WGS) of DNA and gene expression levels from RNA sequencing (RNA-seq) of whole blood in 2622 participants in Framingham Heart Study. We identified 6,778,286 cis-eQTL variant-gene transcript (eGene) pairs at p < 5×10−8 (2,855,111 unique cis-eQTL variants and 15,982 unique eGenes) and 1,469,754 trans-eQTL variant-eGene pairs at p < 1e-12 (526,056 unique trans-eQTL variants and 7,233 unique eGenes). In addition, 442,379 cis-eQTL variants were associated with expression of 1518 long non-protein coding RNAs (lncRNAs). Gene Ontology (GO) analyses revealed that the top GO terms for cis-eGenes are enriched for immune functions (FDR < 0.05). The cis-eQTL variants are enriched for SNPs reported to be associated with 815 traits in prior GWAS, including cardiovascular disease risk factors. As proof of concept, we used this eQTL resource in conjunction with genetic variants from public GWAS databases in causal inference testing (e.g., COVID-19 severity). After Bonferroni correction, Mendelian randomization analyses identified putative causal associations of 60 eGenes with systolic blood pressure, 13 genes with coronary artery disease, and seven genes with COVID-19 severity. This study created a comprehensive eQTL resource via BioData Catalyst that will be made available to the scientific community. This will advance understanding of the genetic architecture of gene expression underlying a wide range of diseases.


Introduction
Over the past decade, genome-wide association studies (GWAS) have revolutionized understanding of the genetic architecture of complex traits. 1 To date, GWAS have reported more than 59,000 associations (at p < 5×10 − 8 ) between common genetic variants and numerous phenotypes (GWAS Catalog, v1.0.2). 2 Yet, despite the clear success of GWAS, most single-nucleotide polymorphisms (SNPs) identi ed in GWAS reside in non-coding regions [3][4][5] and do not illuminate causal mechanisms underlying SNP-trait associations. 5 We posit that many of these trait-associated non-coding SNPs are likely to be involved in the regulation of gene expression.
Expression quantitative trait locus (eQTL) analysis seeks to identify genetic variants that affect the expression of local (cis) or distant (trans) genes (eGenes).
Until recently, eQTL analysis has relied on high throughput microarray technologies and spawned a wave of genome-wide eQTL studies 6-11 including a recent study from our group. 12 These studies aided the understanding of the functional relevance of many GWAS results. Importantly, a hypothesis-free genome-wide eQTL approach permits the identi cation of new putatively functional loci without requiring previous knowledge of speci c regulatory regions.
Most previous eQTL analyses were limited by small sample sizes and by the imprecision of microarrays. Newer technologies of RNA sequencing (RNA-seq) and whole genome sequencing (WGS) of DNA add greater precision and relevance to eQTL analyses. In conjunction with the National Heart, Lung, and Blood Institute's (NHLBI) Trans-Omics for Precision Medicine (TOPMed) Program, 13 the Framingham Heart Study (FHS) has obtained whole genome sequencing (WGS) in ~ 6100 study participants to help understand the molecular basis of heart, lung, blood, and sleep disorders and to advance precision medicine. Among FHS participants with WGS, RNA-seq was obtained in 2622 participants. We conducted genome-wide eQTL analyses using high-precision genotypes obtained via WGS and gene expression levels from RNA-seq of whole blood. The primary objectives of this study were three-fold. Firstly, it sought to provide a scienti c resource of cis and trans gene-level eQTL data to facilitate understanding of the genetic architecture of gene expression traits. Secondly, it was aimed to provide eQTL data for long noncoding RNAs (lncRNAs) that were not captured in prior array-based eQTL studies. Thirdly, it attempted to demonstrate the utility of the eQTL resource in causal inference analyses.

Gene Ontology analyses
We identi ed 100 signi cant GO terms for the top 1000 cis-eGenes at FDR < 0.05. Of these Go terms, there were 58 for Biological Process, 31 for Cellular Component, and 11 for Molecular Function (Supplemental Table 5). Of note, the top GO terms appeared to be related to immune functions. For example, the top two Biological Processes are "leukocyte degranulation" (FDR = 1e-6) and "myeloid leukocyte mediated immunity" (FDR = 2e-6) and the top two Cellular Components are cytoplasm (FDR = 3e-6) and MHC protein complex (FDR = 6e-6). The top 1000 top trans-eGenes gave rise to 75 signi cant (FDR < 0.05) GO terms including 37 for Biological Process, 32 for Cellular Component, and 6 for Molecular Function. The top GO terms for the top 1000 trans-eGenes were enriched in pathways and molecular functions related to immune functions (Supplemental Table 5).

Mendelian randomization analysis
We performed two-sample MR to test for potential causal association of the cis-eGenes with SBP, CAD, and COVID-19 severity. We found 1558 genes containing at least one eQTL variant (median 29; interquartile range [IQR] 6, 88) that coincided with variants from GWAS of SBP (p < 5e-8). 16 After Bonferroni correction for multiple testing, MR identi ed putative causal associations for 60 genes with SBP (i.e., p < 0.05/1558) ( For CAD, 173 genes contained at least one eQTL variant [median 5; IQR (2, 18) that also were associated with CAD in GWAS. 17 Thirteen genes showed putative causal associations with CAD (i.e., p < 0.05/173) ( rs10774671, at exon 7 of OAS1 for which the "G" allele leads to a "prenylated" protein that is protective against severe COVID. 20 Additional MR analysis using rs10774671 as the instrumental variable demonstrated that splice variation of OAS1 is also causal for COVID-19 severity (p = 4e-6).

Replication analyses
Of the reported 10,914 cis-eQTL-eGene pairs from the study by Battle et al. 21 (FDR < 0.05), 21 6782 (62%) pairs displayed p < 5e-8 in the present study. The average proportion of variance explained by these 6782 cis-eQTL variants in respective genes was 0.11 (Supplemental Table 10). Of the 269 trans-eQTL-eGene pairs (FDR < 0.05) reported by Battle et al. 21 47 (18%) pairs displayed p < 1e-12 in the current study. The average proportion of variance explained by these 47 trans-eQTL variants in respective genes was 0.076. Of note, all 47 trans-eQTL variants and respective trans-eGenes are located on the same chromosomes (Supplemental Table 11). The average distance between these trans-eQTL variants and respective trans-eGenes is within 22 Mb.
We conducted additional replication analysis for the cis-eQTL variant-eGene pairs generated from 8,372,247 SNPs and 20,188 gene transcripts that were common to our study (n = 2622 participants) and to GTEx(6) (n = 755 participants) (Supplemental Fig. 1). At p < 5e-8, we identi ed 1,080,485 cis-eQTL variant-eGene pairs in GTEx and 3,852,182 pairs in our study; of these, 951,085 pairs (88% of pairs in GTEx) displayed the same effect direction as in our larger study.

Discussion
We leveraged WGS and RNA-seq data from 2,622 FHS participants to create a powerful scienti c resource of eQTLs. We identi ed signi cant unique cis-eQTL variants-eGene pairs (n = 2,855,111 unique variants with cis-15,982 eGenes) and 526,056 unique trans-eQTL variants-eGene pairs (526,056 unique variants and unique 7,233 trans-eGenes. A large proportion of reported cis-eQTL variant-eGene pairs were replicated with directionally concordant in our study including 88% of cis-variant-eGene pairs from GTEx. Consistent with our previous study and others, 7-12,22,23 90% of eQTL variants identi ed in the present study are located in within 1 Mb of the corresponding cis-eGene and 83% are within 100 kb of the TSSs of the corresponding eGene. While the majority of (85% of cisand 96% of trans-) lead eQTL variants explained only a small proportion (R 2 < 0.2) of interindividual variation in expression of the corresponding eGenes, 15% of lead cis-eQTL variants and 4% of lead trans variant explained 20% or more of interindividual variation in expression of the corresponding eGenes 24 . Additionally, eQTL variants were enriched (p < 0.0001) in disease-associated SNPs identi ed by GWAS. We further demonstrated the utility of our eQTL resource for conducting causal inference testing.
Our MR analyses revealed putatively causal relations of gene expression to several disease phenotypes including SBP, CAD, and COVID-19 severity. Taken together, the comprehensive eQTL resource we provide can advance understanding of the genetic architecture of gene expression underlying a wide variety of diseases. The interactive and browsable eQTL resource will be posted to the National Heart, Lung, and Blood Institute's BioData Catalyst site and will be freely accessible to the scienti c community.
Our study expands current knowledge by creating an accessible and browsable resource of eQTLs based on WGS and RNA-seq technologies. It also includes eQTLs for lncRNAs that were not reported in prior eQTL studies that used array-based expression pro ling. Over the past decade, accumulating evidence shows that lncRNAs are widely expressed and have key roles in gene regulation. 25,26 It is estimated that the human genome contains 16,000 to 100,000 lncRNAs. 25 We identi ed 447,598 cis-eQTL variants for 1518 cis-lncRNAs and 121,241 trans-eQTLs for 475 trans-lncRNAs (Supplemental Tables 3 &4). In addition, we identi ed six lncRNAs that showed putative causal associations with SBP. However, the functions of these six lncRNAs remain to be determined.
Thus, our novel eQTL database may also help in the study of non-protein-coding RNAs in relation to health and disease.
As a proof of concept of the application of the eQTL resource, we performed MR analyses on a small number of cardiovascular traits and COVID-19 severity and demonstrated that the eQTL database can identify promising candidate genes with evidence of putatively causal relations to disease that may merit functional studies. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread across the globe and caused millions of deaths since it emerged in 2019. Recent GWAS of COVID-19 susceptibility and severity [27][28][29] have identi ed SNPs in several loci on chromosomes 3, 9 and 21. 30 Using our eQTL resource in conjunction with COVID-19 GWAS, we conducted MR analyses that identi ed seven genes, including OAS1 and IFNAR2, as putatively causal for COVID-19 severity. The OAS1/2/3 cluster has been identi ed as a risk locus for COVID-19 severity. 27 . This area harbors a protective haplotype of approximately 75 kilo-bases (kb) at 12q24.13 among individuals of European ancestry. 19 A recent study identi ed an alternative splicing variant, rs10774671, at exon 7 of OAS1 for which the protective allele "G" leads to a more active OAS1 enzyme. 20 Our MR results suggest that both the OAS1 gene expression level and its splice variation are causal for COVID-19 severity.
The IFNAR2 gene encodes a protein in the type II cytokine receptor family. Mutations in IFNAR2 are associated with Immunode ciency and measles virus susceptibility and play an essential and a narrow role in human antiviral immunity. 31 A recent study further showed that loss-of-function mutations in IFNAR2 are associated with severe COVID-19. 32 These studies, considered alongside our MR results provide evidence of a causal role of IFNAR2 expression in severe COVID-19 infection.
This study has several noteworthy limitations. This study included White participants of European ancestry who were middle-aged and older; therefore, the eQTLs identi ed may not be generalizable to other races or age ranges. The current RNA-seq platform included ~ 7700 lncRNAs, which is a modest subset of all lncRNAs in the human genome. 25 We used MR analyses to infer causal relation of genes to disease traits. MR analysis is predicated on a set of critical assumptions that may not be testable in the setting of eQTL analysis. 33,34 Replication of our eQTL ndings is warranted in studies with larger sample sizes and more diverse populations.
Our study also has several strengths. The advent of high-throughput RNA sequencing technology provides an unparalleled opportunity to accelerate understanding of the genetic architecture of gene expression. Our study extends and expands the existing literature by identifying novel eQTLs based on WGS and RNA-seq. We demonstrate the potential applications of a vast eQTL resource by analyzing the concordance of eQTL variants with SNPs from GWAS of several disease phenotypes followed by causal inference analyses that identi ed promising disease-related genes that may merit functional studies.
We created an open and freely accessible eQTL repository that can serve as a promising scienti c resource to better understand of the genetic architecture of gene expression and its relations to a wide variety of diseases.

Study participants
This study included participants from the FHS Offspring [10] and Third Generation cohorts [11]. Blood samples for RNA seq were collected from Offspring participants who attended the ninth examination cycle (2011-2014) and the Third Generation participants who attended the second examination cycle (2008)(2009)(2010)(2011). Protocols for participant examinations and collection of genetic materials were approved by the Institutional Review Board at Boston Medical Center. All participants provided written, informed consent for genetic studies. All research was performed in accordance with relevant guidelines/regulations. Isolation of RNA from whole blood and RNA-seq Peripheral whole blood samples (2.5 mL) were collected from FHS participants (Offspring participants at the ninth examination cycle and the Third Generation participants at the second examination cycle) using PAXgene™ tubes (PreAnalytiX, Hombrechtikon, Switzerland), incubated at room temperature for 4 hours for RNA stabilization, and then stored at − 80°C until use. Total RNA was isolated using a standard protocol using a PAXgene Blood RNA Kit at the FHS Genetics Laboratory (FHS Third Generation cohort) and the TOPMed contract laboratory at Northwest Genomics Center (Offspring cohort). Tubes were allowed to thaw for 16 hours at room temperature. White blood cell pellets were collected after centrifugation and washing. Cell pellets were lysed in guanidinium-containing buffer. The extracted RNA was tested for its quality by determining absorbance readings at 260 and 280 nm using a NanoDrop ND-1000 UV spectrophotometer. The Agilent Bioanalyzer 2100 micro uidic electrophoresis (Nano Assay and the Caliper LabChip system) was used to determine the integrity of total RNA.
All RNA samples were sequenced by an NHLBI TOPMed program 13 reference laboratory (Northwest Genomics Center) following the TOPMed RNA-seq protocol. All RNS-seq data were processed by University of Washington. The raw reads (in FASTQ les) were aligned using the GRCh38 reference build to generate BAM les. RNA-SeQC 35 was used for processing of RNA-seq data by the TOPMed RNA-seq pipeline to derive standard quality control metrics from aligned reads. Gene-level expression quanti cation was provided as read counts and transcripts per million (TPM). GENCODE 30 annotation was used for annotating gene-level expression.

Whole blood cell counts
Whole blood cell counts include white blood cell (WBC) count, red blood cell count, platelet count, and WBC differential percentages (neutrophil percent, lymphocyte percent, monocyte percent, eosinophil percent, and basophil percent). Contemporaneously measured blood cell counts were available in 2094 (80%) of the 2622 FHS participants used in eQTL analyses. We performed partial least squares (PLS) prediction method 36 with three-fold cross-validation (2/3 samples for training and 1/3 for validation) to impute these blood cell components using gene expression from RNA-seq. Prediction accuracy (R-squared) varied across blood component: WBC: 0.58, platelet: 27%, neutrophil percentage: 82%, lymphocyte percentage: 85%, monocyte percentage: 77%, eosinophil percentage: 87%, basophil percentage: 32%. Because 80% of the participants in this study had directly measured cell count variables and only 20% received imputed variables, we used the measured (in 2094 participants) and predicted (in 528 participants) blood cell components as covariates in regression models for eQTL analyses.
RNA-seq quality control, and data adjustment To minimize confounding, expression residuals were generated by regressing transcript expression level on age, sex, measured or predicted blood cell count and differential cell proportions, and genetic principal components. Principal component (PC) analysis is a technique for reducing the dimensionality in large data sets. 37 It has been widely used in regression analyses to minimize unknown confounding. We included ve PCs computed from FHS genotype pro les to account for population strati cation. We also included 15 PCs computed from the transcriptome pro le to account for unknown confounders that may affect gene expression. In addition, we adjusted for a relatedness matrix, and technical covariates including year of blood collection, batch (sequencing machine and time, plate and well), and RNA concentration.

Whole genome sequencing
Whole genome sequencing of genomic DNA from whole blood was conducted in ~ 6,000 FHS participants as part of NHLBI's TOPMed program. 13 Standard procedures were used to obtain DNA fragmentation and library construction. Sequencing was performed by a TOPMed reference laboratory (the Broad Institute of MIT and Harvard) using Hi Seq X with sequencing software HiSeq Control Software (HCS) version 3.3.76, then analyzed using RTA2 (Real Time Analysis). The DNA sequence reads were aligned to a human genome build GRCH38 using a common pipeline across all TOPMed WGS centers. A sample's sequence was considered complete when the mean coverage of nDNA was ≥ 30x. This analysis used genetic variants generated from TOPMed Freeze 10a. 13 Association analyses of expression levels with SNPs We performed association analyses of expression levels with genome-wide SNPs with minor allele frequencies (MAFs) ≥ 0.01. In a simple regression model, a SNP was used as an independent variable and the residuals of a transcript expression level was used as the dependent variable. All analyses were performed on the NIH-supported STRIDES cloud infrastructure. A graphical Processing Unit (GPU)-based program 12 was used to facilitate computation. Effect sizes, standard error, partial R-squared, and p-values for all SNP-gene expression pairs were stored to enable complete lookups and to facilitate later meta-analysis.
In this study, we de ned cis-eQTLs as targeting genes within 1 Mb of their transcription start site (TSS). Trans-eQTLs referred to those that were beyond of 1 Mb of the TSSs of the eGenes on the same chromosome or those on the different chromosomes of the eGenes. A signi cant cis-eQTL of an eGene was identi ed if a SNP within 1 Mb of that gene was associated with expression of a transcript of that gene at P < 5x10 − 8 . A signi cant trans-eQTL was de ned as a SNP beyond 1 Mb that gave rise to P < 1x10 − 12 in association a gene.

Gene Ontology analyses
We selected the single, most signi cant eQTL variant (i.e. lead variant) for each eGene (for the gene level analysis) from cisand trans-eQTL results separately. The eGenes annotated to the selected lead cis and trans eQTL variants were matched into Entrez IDs. We used the "goana" function from the "limma" package 38 to test for over-representation of gene ontology (GO) terms or KEGG pathways applied to the top 1000 eGenes. We used FDR < 0.05 to report GO terms including Biological Process, Cellular Component, and Molecular Function.

Enrichment analyses using GWAS Catalog
We linked the eQTL variants with SNPs from the GWAS Catalog 2 (data downloaded on October 22, 2021), which included 243,618 entries for 2,960 mapped traits at p < 5e-8. Cisand trans-eQTL variants were analyzed separately. Unique SNP RS IDs were used for enrichment analysis with Fisher's test. FDR < 0.05 was used for signi cance.
Correlation analysis of selected lncRNA and protein coding genes For lncRNAs that were in the top 25 cis-eQTL variant-eGene pairs, we performed partial Pearson correlation analyses between the expression level of the lncRNA and its nearby protein coding gene, adjusting for the same set of covariates that were included in eQTL analysis. We performed random sampling of 1000 genes 500 times to derive null distributions of partial Pearson correlation of these gene pairs. We calculated an empirical p-value to evaluate whether the partial Pearson correlation coe cient between the expression level of an lncRNA and its nearby protein coding gene was signi cantly higher than the average partial Pearson correlation coe cient from randomly selected gene pairs. The empirical p-value was calculated as the proportion of partial Pearson correlation coe cients that were more extreme than the correlation coe cient of an lncRNA and its nearby protein coding gene.

Mendelian randomization analysis
We conducted Mendelian randomization (MR) to demonstrate the application of the eQTL resource in causal inference analysis. We tested for potential causal association of the cis-eGenes with SBP, coronary artery disease (CAD), and COVID-19 severity. SBP-associated SNPs were obtained from GWAS of over  27 We performed two-sample MR analyses 34 using the TwoSampleMR R package. 39 The instrumental variables (IVs) were independent cis-eQTL variants (LD r 2 < 0.1) from this study. The primary analysis used the inverse variance weighted (IVW) method. We also assessed heterogeneity of the IVs in each gene and conducted sensitivity analysis using the MR-Egger method to test for potential horizontal pleiotropy. We also performed the median-based method 40 and mode-based method 41 when heterogeneity was present in MR analyses due to outliers among the IVs 42 . We reported putative causal genes if Bonferroni correction p < 0.05/n (n is the number of genes tested).

Replication analyses
A previous study reported 10,914 cis-eQTL variant-eGene pairs and 269 trans pairs (FDR < 0.05) through RNA-sequencing of 922 individuals. 21 We performed replication analyses using the reported cis-and trans-eQTL variant-eGene pairs in conjunction with the pairs in the present study. 21 We also used the cis-eQTL database generated from GTEx whole blood (version 8) (https://www.gtexportal.org/home/datasets) for replication of our cis-QTL ndings. Whole genome sequencing and RNA-seq were conducted in whole blood of 755 samples in GTEx. The replication was only performed using the cis-eQTL-variant-eGene pairs generated by 8,372,247 SNPs and 20,188 gene transcripts that were found in common between our study and GTEx. Because this study was aimed to provide eQTL resource for the broad scienti c community, we present replication results using both p < 5e-8 and p < 1e-4 for replicating cis-eQTL variant-eGene pairs.

Declarations
Disclaimer: The views and opinions expressed in this manuscript are those of the authors and do not necessarily represent the views of the National Heart, Lung, and Blood Institute, the National Institutes of Health, or the U.S. Department of Health and Human Services.
Data availability and materials: The datasets analyzed in the present study are available at the dbGAP repository phs000007.v32.p13. The datasets analyzed in the present study are available at the dbGAP repository phs000007.v32.p13 (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi? study_id=phs000007.v30.p11). Variance in eGenes explained by signi cant cis-eQTLs in relation to the distance of signi cant cis-eQTLs to the transcription start site of the cis-gene.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.