Pan-ancestry exome-wide association analyses of COVID-19 outcomes in 586,157 individuals.

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) causes coronavirus disease 2019 (COVID-19), a respiratory illness that can result in hospitalization or death. We used exome sequence data to investigate associations between rare genetic variants and seven COVID-19 outcomes in 586,157 individuals, including 20,952 with COVID-19. After accounting for multiple testing, we did not identify any clear associations with rare variants either exome wide or when specifically focusing on (1) 13 interferon pathway genes in which rare deleterious variants have been reported in individuals with severe COVID-19, (2) 281 genes located in susceptibility loci identified by the COVID-19 Host Genetics Initiative, or (3) 32 additional genes of immunologic relevance and/or therapeutic potential. Our analyses indicate there are no significant associations with rare protein-coding variants with detectable effect sizes at our current sample sizes. Analyses will be updated as additional data become available, and results are publicly available through the Regeneron Genetics Center COVID-19 Results Browser.

Main text

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) causes coronavirus disease 2019 (COVID-19). COVID-19 ranges in clinical presentation from asymptomatic infection to flu-like illness with respiratory failure, hyperactive immune responses, and death.3, 4, 5 Known risk factors for severe disease include male sex, older age, ancestry, obesity, and underlying cardiovascular, renal, and respiratory diseases,6, 7, 8, 9 among others. Since the start of the SARS-CoV-2 pandemic, host genetic analysis of common genetic variation among SARS-CoV-2 patients have identified at least 15 genome-wide significant loci that modulate COVID-19 susceptibility, including variants in/near LZTFL1, IFNAR2, and DPP9.10, 11, 12, 13, 14 However, to date, there has been no exome-wide assessment of the contribution of rare coding genetic variation to COVID-19 disease susceptibility or severity through large population-based exome-wide association analyses.

To identify rare variants (RVs, minor allele frequency [MAF] < 1%) associated with COVID-19 susceptibility and severity, we received approval from institutional review boards (supplemental methods) and analyzed exome-wide sequencing data for 586,157 consented individuals from three studies (Geisinger Health System [GHS], Penn Medicine BioBank [PMBB], and UK Biobank [UKB]) across five continental ancestries (African, Admixed American, European, East Asian, and South Asian; Table S1). Of these, 20,952 had COVID-19, and among those, 4,928 (23.5%) were hospitalized and 1,304 (6.2%) had severe disease (i.e., requiring ventilation or resulting in death; Table S2). Using these data, we tested the association between RVs and seven COVID-19 outcomes: five related to disease susceptibility and two related to disease severity among individuals with COVID-19 (Table S3). In a separate paper, we used these same phenotypes to validate the association with common risk variants reported in previous COVID-19 genome-wide association studies (GWASs),10, 11, 12 , thus demonstrating that our phenotypes are calibrated with those used in other studies.

For each phenotype, exome-wide association analyses were performed separately in each study and ancestry via REGENIE, testing individual RVs (∼7 million) and a burden of RVs in 18,886 protein-coding genes. The genomic inflation factor (λGC) for RVs was often <1 in individual studies, caused by a large proportion of variants having a minor allele count (MAC) of 0 in affected individuals (Table S4). In meta-analyses across studies and ancestries, we found no RV associations at a conservative p < 9.6E−10, which corresponds to a Bonferroni correction for the number of variants and traits tested. At a less conservative significance threshold of p < 5E−8, we found eight genes with RV associations (Table 1 ), of which, we highlight two with an established role in anti-viral responses. First, we highlight an association between higher risk of severe COVID-19 and a burden of ultra-rare (MAF < 0.001%) predicted loss-of-function (pLoF) and missense variants in the toll-like receptor 7 gene (TLR7; p = 4E−8; OR = 4.53; 95% CI = 2.64–7.77), consistent with relatively small exome-sequencing studies of males with severe COVID-19. , TLR7 encodes a single-stranded viral RNA sensor that recognizes coronaviruses, including SARS-CoV-1, MERS, and most likely SARS-CoV-2, and that activates the type-1 interferon pathway in COVID-19. Second, we highlight an association between higher risk of COVID-19 and an ultra-rare missense variant in ZC3HAV1 (rs769102632:A, MAF = 0.002%; p = 3E−8; OR = 26.7; 95% CI 8.37–85.38; Figure S1), a gene that encodes a zinc finger antiviral protein , that inhibits SARS-CoV-2 replication, potentially by upregulating type I interferon responses. Given the potential significance of this finding, we attempted to replicate the ZC3HAV1 rs769102632:A association in an additional 6,223 individuals with COVID-19 with exome or whole-genome sequence data generated as part of the GenOMICC (n = 4,851), Columbia University COVID-19 Biobank (n = 1,152), and Biobanque Quebec (n = 220) studies. We found no carriers for this variant in these additional COVID-19 cases (Table S5) when we expected about four given the observed allele frequency in cases in our study (three and one carriers expected in individuals of African and European ancestry, respectively). Given these findings, we conclude that it is unlikely that there is a true association between rs532051930 and COVID-19 risk. Similarly, the association with a promoter variant in EEF2 that we reported in an earlier version of these analyses was considerably attenuated (from p = 6E−9 to 3E−6), consistent with a false-positive association.

Table 1

Top associations between COVID-19 outcomes and protein-coding rare variants (p < 5E−8)

Gene	Varianta	Variant effect	Odds ratio (95% CI)	p value	N affected individuals with 0|1|2 copies of effect allele	N control individuals with 0|1|2 copies of effect allele	Effect allele frequency	Heterogeneity p value
COVID-19 positive versus COVID-19 negative or unknown
ZC3HAV1	rs769102632	missense	26.72 (8.37, 85.38)	2.95E−8	13,950|7|0	401,218|8|0	0.00002	0.9517
FLNB	rs1256764500	missense	26.6 (8.25, 85.77)	3.97E−8	18,616|7|0	500,616|8|0	0.00001	0.4354
COVID-19 positive versus COVID-19 negative
DISP3	burden	pLoF and deleterious missense with MAF < 10−3	1.88 (1.51, 2.34)	2.26E−8	20,727|145|0	74,172|301|0	0.00234	0.9972
COVID-19 hospitalized versus COVID-19 negative or unknown
WDR78	rs754119466	splice region	49.21 (13.61, 177.85)	2.81E−9	3,619|6|0	392,658|24|0	0.00004	1
TES	rs761377603	missense	38.91 (10.75, 140.9)	2.44E−8	4,555|5|0	511,328|23|0	0.00003	0.6601
MARK1	burden	pLoF variants with MAC = 1	40.19 (10.9, 148.1)	2.86E−8	4,473|5|0	530,595|34|0	0.00004	0.4035
SHC2	rs2287960	stop gained	42.94 (11.17, 165.02)	4.42E−8	4,237|5|0	483,826|17|0	0.00002	0.6742
COVID-19 severe versus COVID-19 negative or unknown
TLR7b	burden	pLoF and missense variants with MAF < 10−5	4.53 (2.64, 7.77)	4.28E−8	1,266|1|7	517,523|383|123	0.00062	0.7188

MAF, minor allele frequency; MAC, minor allele count; CI, confidence interval.

Effect allele for individual variants was rs769102632:A, rs1256764500:G, rs754119466:G, rs761377603:T, and rs2287960:T. For burden tests, individuals were considered to have 0 copies of the effect allele if they were homozygous for the reference allele for all variants included in the burden test, 1 copy of the effect allele if they were heterozygous for at least one variant, and 2 copies if they were homozygous for the alternate allele for at least one variant.

TLR7 is located on the X chromosome. Hemizygous males are included in the N of individuals with two copies of the effect allele.

Next, we addressed the possibility that associations with protein-coding RVs might help pinpoint target genes of common risk variants identified in GWASs of COVID-19. To this end, we focused on 281 genes located within 500 kb of the 15 common risk variants identified by the COVID-19 Host Genetics Initiative (HGI) and asked whether there was any evidence for association between our five COVID-19 susceptibility outcomes and a burden of RVs in any of these genes. We considered associations with pLoF variants alone (M1 burden test) or pLoF together with deleterious missense variants (M3 burden test). No associations surpassed the Bonferroni significance threshold of 3.5E−6, which accounts for the 14,050 gene burden tests performed (281 genes × two burden tests × five allele frequency cut-offs × five susceptibility phenotypes; Table S6). As such, at current sample sizes, RV associations do not point to potential effector genes underlying associations between common variants and COVID-19.

We then examined the association with 13 genes in the interferon pathway, given a previous report that deleterious RVs in these genes may be implicated in severe clinical outcomes. Specifically, we examined whether there was any evidence for association between the COVID-19 hospitalization phenotype (4,928 affected individuals versus 558,763 control individuals) and the burden of rare (MAF < 0.1%, as reported by Zhang et al.) pLoF variants (M1 burden test) or pLoF plus deleterious missense variants (M3 burden test) in these 13 genes. There were no significant associations with any gene, either individually or on aggregate (all burden tests with p > 0.05; Table 2 ). Further, these results were unchanged when testing severe cases of COVID-19 (n = 1,304) or when restricting the burden tests to include variants with an MAF < 1% or singleton variants (Table S7). Therefore, in alignment with a similar report, we also found no evidence for an association between RVs in these 13 interferon-signaling genes.

Table 2

Burden associations among interferon signaling genes

Variants included in burden test	Gene	Odds ratio (95% CI)	p value	N affected individuals with RR|RA|AA genotypea	N control individuals with RR|RA|AA genotypea	AAF	Heterogeneity p value
pLoF, MAF < 0.1%	IFNAR1	1.46 (0.51, 4.17)	0.4786	4,775|5|0	549,164|374|0	0.00034	0.9111
	IFNAR2	1.96 (0.91, 4.19)	0.0844	4,920|8|0	558,068|695|0	0.00062	0.0964
	IKBKGb	0.51 (0.04, 6.57)	0.6048	4,394|0|0	500,582|32|10	0.00005	0.9584
	IRF3	0.91 (0.39, 2.11)	0.8293	4,924|3|1	558,279|483|1	0.00043	0.6339
	IRF7	1.15 (0.57, 2.31)	0.6975	4,920|8|0	557,892|871|0	0.00078	0.5267
	IRF9	0.36 (0.02, 6.96)	0.5024	4,478|0|0	530,571|58|0	0.00005	0.9996
	STAT1	0.36 (0.01, 19.89)	0.6207	4,394|0|0	500,584|40|0	0.00004	0.9996
	STAT2	0.36 (0.07, 1.91)	0.2311	4,644|0|0	541,214|144|0	0.00013	1.0000
	TBK1	0.36 (0.04, 3.13)	0.3553	4,478|0|0	530,539|90|0	0.00008	0.9995
	TICAM1	0.81 (0.14, 4.73)	0.8160	4,477|1|0	530,454|175|0	0.00016	0.7587
	TLR3	1.56 (0.47, 5.13)	0.4656	4,924|4|0	558,457|306|0	0.00027	0.7039
	TRAF3	0.37 (0.0, 217.91)	0.7576	4,394|0|0	500,597|27|0	0.00003	1.0000
	UNC93B1	0.77 (0.28, 2.06)	0.5974	4,641|3|0	540,929|429|0	0.00040	0.9294
	all autosomal genes	0.81 (0.56, 1.18)	0.2709	4,655|23|0	514,810|3,219|0	0.00320	0.9492
pLoF and missense predicted deleterious, MAF < 0.1%	IFNAR1	1.51 (0.71, 3.18)	0.2831	4,918|10|0	557,991|772|0	0.00069	0.8283
	IFNAR2	1.87 (0.88, 3.97)	0.1021	4,920|8|0	558,045|718|0	0.00064	0.0862
	IKBKGb	1.48 (0.18, 12.34)	0.7184	4,393|1|0	500,544|70|10	0.00009	0.6366
	IRF3	0.9 (0.42, 1.92)	0.7778	4,923|4|1	558,128|634|1	0.00057	0.7436
	IRF7	1.15 (0.67, 1.96)	0.6102	4,914|14|0	557,238|1,525|0	0.00137	0.3523
	IRF9	0.36 (0.02, 6.96)	0.5024	4,478|0|0	530,571|58|0	0.00005	0.9996
	STAT1	0.35 (0.08, 1.49)	0.1563	4,762|0|0	547,803|231|0	0.00021	1.0000
	STAT2	1.26 (0.73, 2.2)	0.4089	4,909|19|0	557,153|1,609|1	0.00145	0.7935
	TBK1	1.0 (0.54, 1.85)	0.9951	4,917|11|0	557,567|1,195|1	0.00107	0.6983
	TICAM1	0.8 (0.14, 4.66)	0.8084	4,477|1|0	530,451|178|0	0.00017	0.7558
	TLR3	0.74 (0.49, 1.11)	0.1396	4,911|17|0	556,016|2,745|2	0.00245	0.8319
	TRAF3	1.7 (0.44, 6.62)	0.4431	4,778|2|0	549,284|254|0	0.00023	0.1923
	UNC93B1	0.92 (0.56, 1.5)	0.7309	4,913|15|0	557,079|1,684|0	0.00151	0.9180
	all autosomal genes	0.94 (0.76, 1.17)	0.5835	4,590|88|0	507,793|10,233|3	0.00990	0.5285

Association between the phenotype COVID-19 positive hospitalized versus COVID-19 negative or unknown and 13 genes (12 autosomal) related to interferon signaling that were recently reported to contain rare (MAF < 0.1%) deleterious variants in individuals with severe COVID-19. AAF, alternative allele frequency; CI, confidence interval.

RR, individuals who have genotype reference/reference for all variants included in burden test; RA, individuals who have genotype reference/alternate for at least one variant; AA, individuals who have genotype alternate/alternate for at least one variant.

IKBKG is located on the X chromosome. Hemizygous males are included in the N of individuals with two copies of the effect allele.

Lastly, we performed the same analysis for an additional 32 genes that are involved in the etiology of SARS-CoV-2 infection (ACE2, TMPRSS2), encode therapeutic targets for COVID-19 obtained through the ClinicalTrials database (see web resources) (e.g., IL6R, JAK1), or have been implicated in other immune or infectious diseases through GWASs (e.g., IL33). After correcting for 1,600 burden tests performed (32 genes × five traits × five allele frequency thresholds × two burden tests; Bonferroni significance threshold p < 3.1E−5), there were no significant associations with deleterious RVs for this group of therapeutic target genes for COVID-19 (Table S8).

There are caveats to be considered when interpreting results from this study. First, the five continental ancestry groups considered in our analysis included a small number of individuals with admixed ancestry (specifically, those with two continental ancestries with a likelihood > 0.3; see supplemental methods). For example, individuals with admixed African and European ancestry were included in our analysis of African ancestry. This was done to maximize the number and ancestral diversity of the samples included in our analysis and was adequately controlled for in the association analyses carried out with the whole-genome regression approach implemented in REGENIE (test statistics were not inflated). Second, the burden tests we performed were not designed to identify associations with genes that harbor both risk-increasing and risk-lowering rare variants and are expected to provide limited power in these instances. Other approaches have been developed for these situations, such as SKAT/SKAT-O. However, we have not tested the robustness of these alternative burden tests in the context of multi-ancestry meta-analyses, so we opted against applying them in this study. Third, we used a stringent Bonferroni correction to define significance thresholds that account for multiple testing, which are most likely conservative, given the high correlation between traits and burden tests performed.

In summary, we explored the role of rare coding variants on COVID-19 outcomes on the basis of exome-sequence data, capturing genetic variation not assayed by array genotyping or imputation. We did not find any convincing associations with current sample sizes but will continue to expand our analyses and update results periodically at the Regeneron Genetics Center COVID-19 Results Browser (web resources).