COVID-19 clinical outcomes vary immensely, but a patient’s genetic make-up is an important determinant of how they will fare against the virus. While many genetic variants commonly found in the populations were previously found to be contributing to more severe disease by the COVID-19 Host Genetics Initiative, it isn’t clear if more rare variants found in less individuals could also play a role. This is important because genetic variants with the largest impact on COVID-19 severity are expected to be rarely found in the population, and these rare variants require different technologies to be studies (usually whole-exome or whole-genome sequencing). Here, we combined sequencing results from 21 cohorts across 12 countries to perform a rare variant association study. In an analysis comprising 5,085 participants with severe COVID-19 and 571,737 controls, we found that the gene for toll-like receptor 7 (TLR7) on chromosome X was an important determinant of severe COVID-19. Importantly, despite being found on a sex chromosome, this observation was consistent across both sexes.
Host genetics is a key determinant of COVID-19 outcomes. Previously, the COVID-19 Host Genetics Initiative genome-wide association study used common variants to identify multiple loci associated with COVID-19 outcomes. However, variants with the largest impact on COVID-19 outcomes are expected to be rare in the population. Hence, studying rare variants may provide additional insights into disease susceptibility and pathogenesis, thereby informing therapeutics development. Here, we combined whole-exome and whole-genome sequencing from 21 cohorts across 12 countries and performed rare variant exome-wide burden analyses for COVID-19 outcomes. In an analysis of 5,085 severe disease cases and 571,737 controls, we observed that carrying a rare deleterious variant in the SARS-CoV-2 sensor toll-like receptor TLR7 (on chromosome X) was associated with a 5.3-fold increase in severe disease (95% CI: 2.75–10.05, p = 5.41x10-7). This association was consistent across sexes. These results further support TLR7 as a genetic determinant of severe disease and suggest that larger studies on rare variants influencing COVID-19 outcomes could provide additional insights.
IntroductionDespite successful vaccine programs, SARS-CoV-2 is still a major cause of mortality and widespread societal disruption [1,2]. While disease severity has correlated with well established epidemiological and clinical risk factors (e.g., advanced age, obesity, immunosuppression), these do not explain the wide range of COVID-19 presentations [3]. Hence, individuals without one of these known risk factors may have a genetic predisposition to severe COVID-19[4]. These genetic determinants to severe disease can, in turn, inform about the pathophysiology underlying COVID-19 severity and accelerate therapeutics development [5,6].Previous work on COVID-19 host genetics using genome-wide association studies (GWASs) revealed 23 statistically robust genetic loci associated with either COVID-19 severity or susceptibility [7–11]. Given that most GWASs use genetic data obtained from genome-wide genotyping followed by imputation to measure the association between a phenotype and genetic variation, their reliability and statistical power declines as a variant’s frequency decreases, especially at allele frequencies of less than 1%[12]. Ascertainment of rare genetic variation can be improved with sequencing technology [13]. Rare variants are expected to be enriched for larger effect sizes, due to evolutionary pressure on highly deleterious variants, and may therefore provide unique insights into genetic predisposition to COVID-19 severity. Identifying such genes may highlight critical control points in the host response to SARS-CoV-2 infection.Measuring the effect of rare genetic variants on a given phenotype (here COVID-19) is difficult. Specifically, while variants of large effect on COVID-19 are more likely to be rare, the converse is not true, and most rare variants are not expected to strongly impact COVID-19 severity [14]. Therefore, unless large sample sizes and careful statistical adjustments are used, most rare variant genetic associations studies risk being underpowered, and are at higher risk of false or inflated effect estimates if significant associations are found between COVID-19 and genetic loci. This is exemplified by the fact that several rare variant associations reported for COVID-19 have not been replicated in independent cohorts [15–17].Here, we investigated the association of rare genetic variants on the risk of COVID-19 by combining gene burden test results from whole exome and whole genome sequencing. We build off recent work on exome-wide analyses [17] and include close to 5 times the number of severe cases, with a more genetically diverse cohort, to better study the effect of rare variants on COVID-19. To our knowledge, this is the first rare genetic variant burden test meta-analysis ever performed on a worldwide scale, including 21 cohorts, in 12 countries, including all main continental genetic ancestries.ResultsStudy population and outcomeThe final analysis included up to 28,159 individuals infected with SARS-CoV-2, and up to 597,165 controls from 21 cohorts in 12 countries (Fig 1). Most participants were of European genetic ancestry (n = 576,389), but the consortium also included participants of Admixed American (n = 4,529), African (n = 25,465), East Asian (n = 4,716), Middle Eastern (n = 4,977) and South Asian ancestries (n = 9,943). These resulted in a genetically diverse sample of participants (Fig 2). Participating cohorts enrolled patients based on local protocols, and both retrospective and prospective designs were used. Genetic sequencing was also performed locally, and cohorts were provided with a specific framework for quality control analyses, but each were allowed to deviate based on individual needs. Both exome (n = 11 cohorts) and genome sequencing (n = 10 cohorts) were included in the meta-analyses. The mean age of participants was 55.6 years, and 55.9% were females (S3 Table).10.1371/journal.pgen.1010367.g001Fig 1Maps of countries contributing data to the consortium.Sample sizes (cases and controls) for each phenotype were added and represented on the logarithmic scale by each circle. Relative contribution to each phenotype is represented by the three colors. Maps obtained from https://www.naturalearthdata.com/.10.1371/journal.pgen.1010367.g002Fig 2Participant’s genome projection on the first and second genetic principal components of the 1000G reference panel.AFR: African ancestry. AMR: admixed American ancestry. EAS: east Asian ancestry. EUR: European ancestry. MID: middle eastern ancestry. SAS: south Asian ancestry.We studied three separate outcome phenotypes, as previously described by the COVID-19 Host Genetics Initiative (COVID-19 HGI)[8]. Briefly, the outcome cases were defined according to three standard COVID-19 HGI outcomes: A) severe disease: individuals with SARS-CoV-2 infection who died or required invasive respiratory support (extracorporeal membrane oxygenation, intubation with mechanical ventilation, high-flow oxygen support, or new bilevel or continuous positive airway pressure ventilation), B) hospitalisation: individuals with SARS-CoV-2 who died or required hospitalisation, and C) susceptibility to infection: any individual with SARS-CoV-2 infection. These are also referred to as A2, B2, and C2, respectively, in the COVID-19 HGI meta-analyses [8]. For all three phenotypes, controls were all individuals not classified as cases (including population controls with unknown COVID-19 status). The final meta-analyses included up to 5,085 cases and 571,737 controls for the severe disease outcome, 12,304 cases and 590,151 controls for the hospitalisation outcome, and 28,196 cases and 597,165 controls for the susceptibility outcome.Single-variant analysisWe first performed an exome-wide association study using single variants with a MAF (minor allele frequency) higher than 0.1% and an allele count of 6 or more in at least one cohort, with the same additive model and covariates used in the COVID-19 HGI GWAS [8]. Analyses were performed separately by each cohort and each ancestry using Firth regression as applied in the Regenie software [18]. Firth regression is a penalized likelihood regression method that provides unbiased effect estimates even in highly unbalanced case-control analyses [19]. The summary statistics were then meta-analyzed with a fixed effect inverse-variance weighted model within each ancestry, and then with a DerSimonian-Laird random effect model across ancestries.The previously described Neanderthal chromosome 3 locus associated with COVID-19 outcomes [8,20] was also found in all three phenotypes (Figs 3 and S1), with lead variants in the CXCR6 gene for the severe COVID-19 phenotype (rs13059238), and in FYCO1 in the hospitalisation phenotype (rs13069079), and for the LIMD1 gene in the susceptibility phenotype (rs141045534). Reassuringly, each cohort provided summary statistics in the chromosome 3 locus, suggesting that the QC process was working as intended (allowing for sample sizes and number of cases) (S2–S4 Figs).10.1371/journal.pgen.1010367.g003Fig 3Single variant exome-wide association study Manhattan plot (MAF>0.1%).QQ-plot available in the S1 Fig. Black dashed line demarcates the genome-wide significance threshold (p < 5x10-8).Three other loci were found for the hospitalization phenotypes. One at SRRM1 (rs1479489847, OR: 4.17, 95% CI 2.60–6.70, p = 3.25x10-9), IL6R (rs911647797, OR: 6.19, 95% CI: 3.29–11.6, p = 1.45x10-8), and another at cytoskeleton FRMD5 (rs1369031075, OR: 4.06, 95% CI: 2.70–6.11, p = 1.75x10-11). While these loci may hold biological plausibility (especially IL6R, given the use of IL-6 receptor inhibitors in the treatment of COVID-19), these associations were driven by two smaller cohorts (Genentech and Vanda, S5 Fig). However, the SRMM1 locus is located between two stretches of T nucleotides, while both the FRMD5 and the IL-6 loci are within GC rich regions, making variant calling difficult. Hence, these findings will require validation, despite the biological plausibility.Finally, all genetic inflation factors were below 1 (S4 Table). Summary statistics for genome-wide significant variants can be found in S5 Table and QQ-plots can be found in S1 Fig.Burden test definitionGiven the expected paucity of large-effect size rare deleterious variants, strategies have been devised to increase statistical power to test associations between rare variants and biomedically-relevant outcomes. One such strategy is to use burden tests [21], where each variant is collapsed into larger sets of variants, and association is tested between groups of variants and an outcome. Here, we collapsed deleterious variants in each gene and devised the following burden test: for each gene, an individual received a score of 0 if they do not carry any deleterious variant, a score of 1 if they carry at least one non-homozygous deleterious variant, and a score of 2 if they carry at least 1 homozygous deleterious variant. As defined in previous studies on burden testing of rare variants[17,22] deleterious variants were chosen using three masks: 1) “M1" which uses only predicted loss of function variants, 2) “M3” which uses all variants in M1, as well as indels of moderate consequence as predicted by Ensembl [23], and missense variants classified as deleterious in 5 in-silico algorithms (see Methods), and 3) “M4”, which uses all variants in M1 and M3, and also adds all missense variants classified as deleterious in at least 1 of the in-silico algorithms.The analyses were performed separately both for variants with MAF of less than 1%, and for variants of MAF less than 0.1%. We defined MAFs based on a combination of gnomAD [24] MAF annotations, and of cohort-specific common variant exclusion lists. These common variant exclusion lists included variants that achieved a MAF of >1% or >0.1% in at least one study population within the consortium. To reduce the effect of fluctuations due to sampling, a minor allele count (MAC) ≥ 6 in the corresponding study was required for inclusion in the common variant list. Such “blacklists” have been shown to increase statistical power by removing variants at lower risk of being highly deleterious, and it reduces the risk of having cohort-specific false-positive variants being retained on the overall analysis[25]. Hence for each MAF threshold, each cohort removed any variant with a MAF above the threshold in either gnomAD or the corresponding common variant exclusion list.The resulting score (either 0, 1, or 2) for each mask was then regressed on each of our three phenotypes using logistic regression, controlling for age, age^2, sex, sex*age, sex*age2, and 10 common variant (MAF > 1%) genetic principal components (the same covariates as for COVID-19 HGI GWASs [7,8]). Additionally, given that population genetic structure and its confounding effect on phenotypes is different at the rare variant level [26], we also used the first 20 genetic principal components from rare variants (MAF<1%) as covariates in all our analyses. Analyses were otherwise done using the same approach as for single-variant analyses.Exome-wide burden test analyses resultsOur meta-analysis included a total of 18,883 protein-coding genes, and all burden test genetic inflation factors, for all masks, were less than 1 (S4 Table), suggesting that our results were not biased by population stratification and that Firth regression adequately adjusted for unbalanced case-control counts. Using an exome-wide significance p-value threshold of 0.05/20,000 = 2.5x10-6, we found 3 genes associated with one of the COVID-19 phenotypes in at least one mask in our meta-analyses (Table 1 and S6–S12 Figs). Of specific interest, we observed that carrying a predicted loss of function or in-silico highly deleterious missense variant (i.e., mask M3) in the toll-like receptor 7 (TLR7) gene was associated with a 5.3-fold increase (95% CI: 2.7–10.1, p = 5.41x10-7) in odds of severe COVID-19. TLR7 is an important part of the innate viral immunity, encoding a protein that recognizes coronaviruses and other single-stranded RNA viruses, leading to upregulation of the type-1 and type-2 interferon pathway [27]. Results from the severe COVID-19 outcome analyses of TLR7 with other masks also nearly reach our statistical significance threshold, with larger effects found in the M1 mask (OR: 13.6, 95% CI: 4.41–44.3, p = 1.64x10-5) and smaller effect in the M4 mask (OR: 3.12, 95% CI: 1.91–5.10, p = 5.30x10-6), though the latter was balanced by smaller standard errors due to the larger number of cases (3275 cases in M4 vs 1577 in M1), as expected. These findings further support previous reports of TLR7 errors of immunity underlying severe COVID-19 presentations [17,28–31].10.1371/journal.pgen.1010367.t001Table 1Exome-wide significant findings, as well as other TLR7 results (for the severe phenotype only).Note that for Masks M1, all deleterious variants had a MAF<0.1%, and hence both burden tests (MAF<1% and 0.1%) gave the same results. Full results available in S4 Table.GeneMaskPhenotypeMAFBetaStandard ErrorOdds Ratio95% Confidence IntervalP-valueHeterogeneity p-valueN Cases 0|1|2 Burden TestN Controls 0|1|2 Burden TestMeta-Analysis Across AncestriesMARK1M1Severe COVID-19<0.1%3.170.6723.96.5–88.21.89x10-60.8831935|4|0540031|92|0MARK1M1Severe COVID-19<1%3.170.6723.96.5–88.21.89x10-60.8831935|4|0540031|92|0MARK1M1Hospitalisation<0.1%2.510.4812.34.8–31.21.43x10-70.8936132|8|0547943|93|0MARK1M1Hospitalisation<1%2.510.4812.34.8–31.21.43x10-70.8936132|8|0547943|93|0RILPL1M1Severe COVID-19<0.1%3.010.6420.25.8–70.72.42x10-60.9411745|4|0558448|121|0TLR7M3Severe COVID-19<0.1%1.660.335.252.75–10.055.41x10-70.7553101|2|5519047|83|47TLR7M3Severe COVID-19<1%1.630.335.102.67–9.727.48x10-70.7603275|2|5519834|85|47Other TLR7 results for severe phenotypeTLR7M1Severe COVID-19<0.1%2.610.6013.64.14–44.41.64x10-50.8201577|0|2508987|13|11TLR7M1Severe COVID-19<1%2.610.6013.64.14–44.41.64x10-50.8201577|0|2508987|13|11TLR7M4Severe COVID-19<0.1%1.140.253.121.91–5.105.30x10-60.8543275|3|7519616|210|139
TLR7 M4 Severe COVID-19 <1% 1.11 0.24 3.03 1.90–4.85 3.43x10-6 0.956 3273|5|8 521166|221|144
In the meta-analyses, we also found that pLoFs (M1) in MARK1 were associated with a 23.9-fold increase in the odds of severe COVID-19 (95% CI: 6.5–88.2, p = 1.89x10-6), and a 12.3-fold increase in the odds of hospitalisation due to COVID-19 (95% CI: 4.8–31.2, p = 1.43x10-7). While the number of MARK1 pLoFs (M1) found in severe and hospitalized cases was small (MAC = 4 and MAC = 8, respectively), the signal was consistent in our three largest cohorts: UK Biobank, Penn Medicine, and Geisinger Health Services (S9–S12 Figs). MARK1 is a member of the microtubule affinity-regulating kinase family, and is involved in multiple biological processes, chief among which is the promotion of microtubule dynamics [32]. MARK1 has previously been shown to interact with the SARS-CoV-2 ORF9b protein [33], further supporting its potential role in COVID-19. Lastly, our meta-analyses also found marginal evidence for an association between severe COVID-19 and pLoFs (M1) in RILPL1 (OR: 20.2, 95% CI: 5.8–70.7, p = 2.42x10-6), a gene that, like MARK1, is associated with microtubule formation and ciliopathy [34].We then meta-analyzed p-values using the aggregated Cauchy association test [35] (ACAT). ACAT accounts for correlation between test statistics (as is expected here) by treating p-values as Cauchy random variables, and taking their weighted average, which also is Cauchy distributed. With ACAT, the association between TLR7 and severe COVID-19 (p = 1.58x10-6), and between MARK1 and hospitalisation (p = 4.30x10-7) remained exome-significant (Fig 4). Full summary statistics are available in S6 Table.10.1371/journal.pgen.1010367.g004Fig 4Exome burden test ACAT p-value meta-analysis Manhattan plots and QQ plots.QQ-plot available in the S6 Fig. Black dashed line demarcates the Bonferroni significance threshold (p < 0.5/20,000).Finally, we note that for both TLR7 and MARK1, the signal was driven by European ancestry participants. Further, while the larger biobanks contributed to these findings, smaller prospective cohorts also provided cases with rare variants at both genes, highlighting the importance of study design in rare variant association testing (S6–S12 Figs).TLR7 sex stratified analysesGiven that TLR7 is located on the X chromosome, we performed sex-stratified analyses of the severe disease phenotype to determine if the effect was also observed in females. These could only be done for the M3 and M4 masks due to very low number of M1 mask qualifying variants (Fig 5). In both we still see a clear effect among males with a 4.81-fold increase in the odds of severe COVID-19 in M3-variant carriers (95% CI: 2.41–9.59, 5 case carriers, 47 control carriers), and a 3.08-fold increase in M4-variant carriers (95% CI: 1.83–5.20, 7 case carriers, 143 control carriers). In females, we still observed a nominally significant signal in the M3 mask, with a 15.2-fold in odds of severe disease in M3-variant carriers (95% CI: 1.51–153.4). However, an M3-variant was observed in only one female with severe disease (heterozygous) in these analyses (compared to 76 heterozygous controls). In M4 variant, the analyses included 2 female heterozygous carriers (and 203 heterozygous controls), with a 4.86-fold in odds of severe disease (95% CI: 0.43–54.3).10.1371/journal.pgen.1010367.g005Fig 5Sex-stratified TLR7 analyses.Rare variants in interferon-related genes and at previously reported genome-wide significant lociDespite a 7.7-fold increase in number of cases, and a 1,069-fold increase in number of controls, the previously reported associations of genes in the interferon pathway with COVID-19 outcomes [15,16] could not be replicated with either our exome-wide significance threshold (S7 Table) or a more liberal one of p = 0.05/10 = 0.005 (based on Bonferroni correction by the number of genes in the interferon pathway defined in a previous study[15]).We also tested for rare variant associations between GWAS candidate genes from genome-wide significant loci in the COVID-19 HGI GWAS meta-analyses, but observed no exome-wide significant associations (S8 Table). However, at a more liberal Bonferroni threshold of p = 0.05/46 = 0.001 (correcting for the 46 genes in the COVID-19 HGI GWAS associated loci), we observed an increased burden of pLoF (M1) or missense variants (M3 mask) in ABO gene among those susceptible to SARS-CoV-2 infection (Table 2). For example, individuals carrying a pLoF (M1) with MAF<0.1% in ABO were at a 2.34-fold higher risk of having a positive SARS-CoV-2 infection (95% CI: 1.50–3.64, p = 1.6x10-4). The ABO results were driven mainly by European and African ancestry participants (S13 Fig). Note that deleterious variants in ABO often lead to blood groups A and B [36,37], which is consistent with the epidemiological association that non-type-O individuals are at higher risk of COVID-19[38]. However, more work is required to better understand the genetics of this locus as it relates to COVID-19 outcomes. Lastly, missense variants in NSF (mask M4, MAF<1%) were also associated with higher susceptibility to SARS-CoV-2 (OR: 1.48, 95% CI: 1.21–1.82, p = 1.4x10-4), but this association was not present in other masks (S9 Table).10.1371/journal.pgen.1010367.t002Table 2Results of burden tests at genes identified from common variants GWAS in the COVID-19 HGI.Only genes with p<0.05/46 are shown here. Full results available in S8 Table.GeneMaskPhenotypeMAFBetaStandard ErrorOdds Ratio95% Confidence IntervalP-valueHeterogeneity p-valueN Cases 0|1|2 Burden TestN Controls 0|1|2 Burden TestNSFM4Susceptibility<1%0.3950.1041.4841.21–1.821.44x10-40.86625752|127|2585642|1907|5ABOM1Susceptibility<0.1%0.8510.2262.3411.50–3.651.68x10-40.49822778|27|0572310|296|0ABOM1Susceptibility<1%0.7840.2092.191.45–3.301.75x10-40.82623460|34|0574608|364|0ABOM3Susceptibility<1%0.7290.1952.0731.41–3.041.89x10-40.86924455|42|0575051|434|0ABOM1Hospitalisation<0.1%1.330.3953.781.74–8.207.56x10-40.5427859|12|0561642|291|0ABOM1Susceptibility<0.1%0.7360.2222.0881.35–3.239.35x10-40.51223779|29|0572799|320|0Replication in GenOMICCData for the M1 mask for TLR7 and MARK1 in the severe COVID-19 phenotype was then replicated with the GenOMICC cohort [11], a prospective study enrolling critically ill individuals with COVID-19, with controls selected from the 100,000 genomes cohort[39]. Results are shown in Table 3. For TLR7, European ancestry individuals with a pLoF (M1) had a 4.70-fold increase in odds of severe disease (95% CI: 1.58 to 14.0, p = 0.005). In the sample of South Asian ancestry individuals, a pLoF (M1) was associated with a 1.90-fold increase in odds of severe disease, but the 95% confidence interval crossed the null (0.23 to 15.6, p = 0.55), which was likely due to a much smaller sample size than in the European ancestry subgroup (1,202 vs 10,645). Of interest, in both Europeans and South Asians, no pLoFs were observed in either of the control groups.10.1371/journal.pgen.1010367.t003Table 3Replication of M1 mask, severe COVID-19, MARK1 and TLR7 results in the GenOMICC cohort. Note that the same variants were included in both the MAF<1% and MAF<0.1% replication, and the same results were obtained (shown here).GeneAncestryBetaStandard ErrorOdds Ratio95% Confidence IntervalP-valueN Cases 0|1|2 Burden TestN Controls 0|1|2 Burden TestMARK1EUR0.1951.421.210.075–19.70.8915988|1|04655|1|0MARK1SAS1.442.194.210.058–307.60.511787|1|0414|0|0TLR7EUR1.550.5554.701.58–14.00.0055980|1|84566|0|0TLR7SAS0.6401.081.900.230–15.60.552786|0|2414|0|0On the other hand, we could not replicate an effect from MARK1, which demonstrated an OR of 1.21 in European ancestry participants (95% CI 0.075 to 19.7, p = 0.89) and an OR of 4.21 in South Asian ancestry individuals (95% CI 0.058 to 307, p = 0.51).DiscussionWhole genome and whole exome sequencing can provide unique insights into genetic determinants of COVID-19, by uncovering associations between rare genetic variants and COVID-19. Specifically, gene burden tests can be particularly helpful, because they test for coding variants, thereby pointing directly to a causal gene and often suggesting a direction of effect. However, such studies require careful control for population stratification and an adapted analysis method such as burden testing, in order to have enough statistical power to find those associations. In our study, we observed that individuals with rare deleterious variants at TLR7 are at increased risk of severe COVID-19 (up to 13.1-fold increase in odds in those with pLoFs). Although this association was suggested by previous studies [28–30], our study provides the most definitive evidence for the role of TLR7 in COVID-19 pathogenesis, with exome-wide significance for this gene in the discovery phase followed by strong replication in a large independent cohort. TLR7 is a well-studied part of the antiviral immunity cascade and stimulates the interferon pathway after recognizing viral pathogen-associated molecular patterns. Given its location on the X chromosome, it has been hypothesis that it could partly explain the observed COVID-19 outcome differences between sexes [40–42], and to our knowledge, this is the first study to show that even in heterozygous females, this gene can potentially play a role in severe disease. Further, this our results suggest that TLR7 mediated genetic predisposition to severe COVID-19 may be a dominant or co-dominant trait, an observation that cannot be made in cohorts limited to male participants[28,30].We also uncovered a potential role for cellular microtubule disruption in the pathogenesis of COVID-19 and the microtubule network is known to be exploited by other viruses during infections [43]. Indeed, the MARK1 protein has been shown to interact with SARS-CoV-2 in previous in-vitro experiments [33]. Nevertheless, these findings at MARK1 were not replicated in the GenOMICC cohort and will need to be tested in larger cohorts, especially given the small number of highly deleterious variants that we found in our consortium. Lastly, we found single variant associations at IL6R, SRRM1, and FRMD5. While IL6R is is already a therapeutic target [44,45] for COVID-19, and SRRM1 has been reported in a previous pre-print [46], these were found in smaller cohorts and will require replication.To our knowledge, this is the first time a rare variant burden test meta-analysis has been attempted on such a large scale. Our framework allowed for easy and interpretable summary statistics results, while at the same time preventing participant de-identification or any breach of confidentiality that stems from sharing results of rare genetic variant analyses [47]. It also provides important insights into how these endeavours should be planned in the future. First, our burden test operated under the assumption that the effect of any of the deleterious variants on the phenotype would be in the same direction and did not account for compound deleterious variant heterozygosity. This allowed for easier meta-analysis across cohorts, but may have decreased statistical power. Other methods may be needed in future analysis to soften this assumption, though some of these cannot be easily meta-analyzed across multiple cohorts directly from summary statistics (e.g., SKAT-O [48]). Similarly, methods that combine both rare and common variants might also provide additional insights into disease outcomes [31,49]. Second, our results highlight the importance at looking at different categories of variants through different masks to increase sensitivity and specificity of our burden tests. Third, while the largest biobanks contributed the most to the signal observed at TLR7 and MARK1, many of our smaller prospective COVID-19 specific cohorts also contributed to the signal. This further highlights the importance of robust study design to improve statistical power, especially with rare variant associations. Lastly, work remains to be done to standardize sequencing and annotation pipelines to allow comparisons of results easily across studies and cohorts. Here, we provided a pipeline framework to every participating cohort, but there remains room for process harmonization. While the decentralized approach to genetic sequencing, quality control, and analyses allowed for more rapid generation of results, it may come at the cost of larger variance in our estimates. In the future, more sophisticated approaches may be required to increase statistical power of exome-wide rare variant association studies [50].Our study had limitations. First, even if this is one of the world’s largest consortia using sequencing technologies for the study of rare variants, we remain limited by a relatively small sample size. For example, in a recent analyses of UK Biobank exomes, many of the phenotypes for which multiple genes were found using burden tests had a much higher number of cases than in our analyses (e.g. blonde hair colour, with 48,595 cases) [22]. Further, rare variant signals were commonly found in regions enriched in common variants found in GWASs. The fact that ABO and NSF were the only genes from the COVID-19 HGI GWAS that were also identified in our burden test (albeit using a more liberal significance threshold), also suggests a lack of statistical power. Similarly, GenOMICC, a cohort of similar size, was also unable to find rare variant associations using burden tests [11]. However, their analysis methods were different from ours, making further comparisons difficult. Nevertheless, this provides clear guidance that smaller studies looking at the effect of rare variants across the genome are at considerable risk of finding both false positive and false negative associations. Second, many cohorts used population controls, which may have decreased statistical power given that some controls may have been misclassified. However, given that COVID-19 critical illness remains a rare phenomenon [51], our severe disease phenotype results are unlikely to be strongly affected by this. Finally, the use of population control is a long-established strategy in GWAS burden tests [7,8,11,22,52], and the statistical power gain from increasing our sample size is likely to have counter-balanced the misclassification bias.In summary, we reproduced an exome-wide significant association with severe COVID-19 outcomes in carriers of rare deleterious variants at TLR7, for both sexes. Our results also suggest an association between the cellular microtubule network and severe disease, which requires further validation. More importantly, our results underline the fact that future genome-wide studies of rare variants will require considerably larger sample size, but our work provides a roadmap for such collaborative efforts.MethodsEthics statementEach cohort had the following statement to make on ethics:BQC-19Each participant or their legal representative (if the participant was incapable to consent) provided informed consent to the biobank. If a participant regained capacity to give consent, informed consent was obtained again directly from the participant. The study was approved by the Jewish General Hospital and Centre Hospitaler de l’Université de Montréal institutional review boards.Columbia Biobank: Recruitment and sequencing of participants from the Columbia COVID-19 Biobank were approved by the Columbia University Institutional Review Board (IRB) protocol AAAS7370 and the genetic analyses were approved under protocol AAAS7948. A subset of patients was included under a public health crisis IRB waiver of consent specifically for COVID-19 studies if patients were deceased, not able to consent, or if the study team was unable to contact them as per IRB protocol AAAS7370.DeCOIInformed consent was obtained from each participant or the legal representative. DeCOI received ethical approval by the Ethical Review Board (ERB) of the participating hospitals/centres (Technical University Munich, Munich, Germany; Medical Faculty Bonn, Bonn, Germany; Medical Board of the Saarland, Germany; University Duisburg-Essen, Germany; Medical Faculty Duesseldorf, Duesseldorf, Germany)FHoGIDEach participant or their legal representative provided informed consent to the biobank. FHoGID received ethical approval by the Commission cantonale d’éthique de la recherche sur l’être humain.GEN-COVID multicenter study: The patients were informed of this research and agreed to it through the informed consent process. The GEN-COVID is a multicentre academic observational study that was approved by the Internal Review Boards (IRB) of each participating centre (protocol code 16917, dated March 16, 2020 for GEN-COVID at the University Hospital of Siena).GenentechThe protocol was reviewed by the institutional review board or ethics committee at each site. Written informed consent was obtained from all the patients or, if written consent could not be provided, the patient’s legally authorized representative could provide oral consent with appropriate documentation by the investigator. Details on institutional review boards are provided in S9 Table.GenOMICCGenOMICC was approved by the appropriate research ethics committees (Scotland, 15/SS/0110; England, Wales and Northern Ireland, 19/WM/0247). Informed consent was obtained for all participants.Geisinher Health Systems: All subjects consented to participation and the analysis was approved by the Geisinger Institutional Review Board under project number 2006–0258.Helix Exome+ and Healthy Nevada Project COVID-19 Phenotypes: informed consent was obtained for all participants. The Healthy Nevada Project study was reviewed and approved by the University of Nevada, Reno Institutional Review Board (IRB, project 956068–12)Thai Biobank ()Informed consent was obtained for each participant via the biobank. The study was approved by the Institutional Review Board of the Faculty of Medicine, Chulalongkorn University, Bangkok, Thailand (COA No. 691/2021).Japan COVID-19 Task Force: Each participant or their legal representative (if the participant was incapable to consent) provided informed consent to the biobank. Study was approved by the ethical committees of Keio University School of Medicine, Osaka University Graduate School of Medicine, and affiliated institutes.Interval WGSAfter reading study leaflets and participating in a discussion with donor carer staff, eligible donors were asked to complete the trial consent form before giving a blood donation. The National Research Ethics Service (United Kingdom) approved this study.MNM Diagnostics (Polish Covid WGS)All participants, or their guardians/parents for the participants under 18), provided their informed consent before collecting their blood samples. The study was approved by the Institutional Ethics Committee of the Central Clinical Hospital of the Ministry of Interior and Administration in Warsaw, Poland (decision nr: 41/2020 from 03.04.2020 and 125/2020 from 1.07.2020).MSCICThis research protocol was reviewed and approved by the Icahn School of Medicine at Mount Sinai Institutional Review Board (IRB) (STUDY-20-00341). During the height of the SARS-CoV-2 pandemic in New York City, all patients admitted to the Mount Sinai Health System were made aware of the research study by a notice included in their admission paperwork. The notice outlined details of the planned research, potential specimen collection and the opportunity to opt-out of research. Flyers announcing the study were also posted throughout the health system. Given the monumental hurdles of consenting sick and infectious patients in isolation rooms, the IRB allowed for specimen collection to occur prior to obtaining research consent at the time of clinical blood collection. Patients and/or their legally authorized representative provided consent to the research study, including genetic profiling for research and data sharing on an individual level. In a subset of individuals, who were unreachable following hospital discharge, we were unable to obtain written informed consent. In these cases, data cannot be share further. All data used these these analyses were anonymized same as above.Penn medicineRecruitment of PMBB participants was approved under IRB protocol 813913 and supported by Perelman School of Medicine at University of Pennsylvania.POLCOVID-GenomikaAll study participants provided written informed consent and received detailed information on the study and associated risk before enrollment. The study was approved by the Bioethics Committee of the Medical University of Bialystok.Qatar Genome ProgramAll QBB participants signed an Informed Consent Form prior to their participation; QBB study protocol ethical approval was obtained from the Hamad Medical Corporation Ethics Committee in 2011 and continued with QBB Institutional Review Board (IRB) from 2017 onwards and it is renewed on an annual basisSaudi human genome programInformed Consent was provided to each participant or their legal guardian (if the participant could not consent) by the corresponding institute. This study was approved the IRB of each participating hospitals, and the IRB at King Abdullah International Medical Research Centre, Ministry of National Guard–Health Affairs, Riyadh, Ministry of Health, and King Fahad Medical City.Swedish BiobankInformed consent was obtained for all study participants. The study was approved by the National Ethical Review Agency (Sweden) (No. 2020–01623).UK BiobankAll subjects consented to participation. The UK Biobank was approved by the North West Multi-centre Research Ethics Committee (United Kingdom) (11/NW/0382). The work described herein was approved by the UK Biobank under application no. 26041.University of California, Los Angeles biobankEach participant or their legal representative (if the participant was incapable to consent) provided informed consent to the biobank. If a participant regained capacity to give consent, informed consent was obtained again directly from the participant. This study was considered human subjects research exempt because all genetic and electronic health records were de-identified. This study was approved by the UCLA Health Institutional Review Board.Vanda COVID-19All participants consented to WGS. The study was reviewed and approved by Advarra IRB; Pro00043096.COVID-19 outcome phenotypesFor all analyses, we used three case-control definitions: A) Severe COVID-19, where cases were those who died, or required either mechanical ventilation (including extracorporeal membrane oxygenation), high-flow oxygen supplementation, new continuous positive airway pressure ventilation, or new bilevel positive airway pressure ventilation, B) Hospitalized COVID-19, where cases were all those who died or were admitted with COVID-19, and C) Susceptibility to COVID-19, where cases are anyone who tested positive for COVID-19, self-reported an infection to SARS-CoV-2, or had a mention of COVID-19 in their medical record. For all three, controls were individuals who did not match case definitions, including population controls for which case status was unknown (given that most patients are neither admitted with COVID-19, nor develop severe disease [53]). These three analyses are also referred to as analyses A2, B2, and C2 by the COVID-19 Host Genetics Initiative [8], respectively.Cohort inclusion criteria and genetic sequencingAny cohort with access to genetic sequencing data and the associated patient level phenotypes were allowed in this study. Specifically, both whole-genome and whole-exome sequencing was allowed, and there were no limitations in the platform used. There were no minimal number of cases or controls necessary for inclusion. However, the first step of Regenie, which was used to perform all tests (see below), uses a polygenic risk score which implicitly requires that a certain sample size threshold be reached (which depends on the phenotype and the observed genetic variation). Hence, cohorts were included if they were able to perform this step. All cohorts obtained approval from their respective institutional review boards, and informed consent was obtained from all participants. More details on each cohort’s study design and ethics approval can be found in the S3 and S1 Tables.Variant calling and quality controlVariant calling was performed locally by each cohort, with the pre-requisite that variants should not be joint-called separately between cases and controls. Quality control was also performed individually by each cohort according to individual needs. However, a general quality control framework was made available using the Hail software [54]. This included variant normalization and left alignment to a reference genome, removal of samples with call rate less than 97% or mean depth less than 20. Genotypes were set to unknown if they had genotype quality less than 20, depth less than 10, or poor allele balance (more than 0.1 for homozygous reference calls, less than 0.9 for homozygous alternative calls, and either below 0.25 or above 0.75 for heterozygous calls. Finally, variants were removed from if the mean genotype quality was less than 11, mean depth was less than 6, mean call rate less than or equal to 0.8, and Hardy-Weinberg equilibrium p-value less than or equal to 5x10-8 (10−16 for single variant association tests). Details on variant calling and quality control is described for each cohort in the S3 Table.Single variant association testsWe performed single variant association tests using a GWAS additive model framework, with the following covariates: age, age2, sex, age*sex, age2*sex, 10 genetic principal components obtained from common genetic variants (MAF>1%). Each cohort performed their analyses separately for each genetic ancestry, but also restricted their variants to those with MAF>0.1% and MAC>6. Summary statistics were then meta-analyzed using a fixed effect model within each ancestry and using a DerSimonian-Laird random effect model across ancestries with the Metal package [55] and its random effect extension [56]. Lastly, given that multiple technologies were used for sequencing, and that whole-exome sequencing can provide variant calls of worse quality in its off-target regions [57], we used the UKB, GHS, and Penn Medicine whole-exome sequencing variants as our “reference panel” for whole-exome sequencing. Hence, only variants reported in at least one of these biobanks were used in the final single-variant analyses.Variant exclusion listFor the burden tests, we also compiled a list of variants that had a MAF > 1% or > 0.1% in any of the participating cohorts. This list was used to filter out variants that were less likely to have a true deleterious effect on COVID-19, even if they were considered rare in other cohorts, or in reference panels [25]. We created two such variant exclusion lists: one to be used in our burden test with variants of MAF less than 1%, and the other for the analysis with MAF less than 0.1%. In any cohort, if a variant had a minor allele count of 6 or more, and a MAF of more than 1% (or 0.1%), this variant was added to our exclusion list. This list was then shared with all participating cohorts, and all variants contained were removed from our burden tests.Gene burden testsThe following analyses generally followed the methods used by recent literature on large-scale whole-exome sequencing [22] and the COVID-19 HGI [8].The burden tests were performed by pooling variants in three different variant sets (called masks), as described in recent UK Biobank whole-exome sequencing papers by Backman et al.[22] and Kosmicki et al.[17].: “M1” which included loss of functions as defined by high impact variants in the Ensembl database[23] (i.e. transcript ablation, splice acceptor variant, splice donor variant, stop gained, frameshift variant, stop lost, start lost, transcript amplification), “M3” which included all variants in M1 as well as moderate impact indels and any missense variants that was predicted to be deleterious based on all of the in-silico pathogenicity prediction scores used, and “M4” which included all variants in M3 as well as all missense variants that were predicted to be deleterious in at least one of the in-silico pathogenicity prediction scores used. For in-silico prediction, we used the following five tools: SIFT [58], LRT [59], MutationTaster[60], PolyPhen2[61] with the HDIV database, and PolyPhen2 with the HVAR database. Protein coding variants were collapsed on canonical gene transcripts.Once variants were collapsed into genes in each participant, for each mask, genes were given a score of 0 if the participant had no variants in the mask, a score of 1 if the participant had one or more heterozygous variant in this mask, and a score of 2 if the participant had one or more homozygous variant in this mask. These scores were used as regressors in logistic regression models for the three COVID-19 outcomes above. These regressions were also adjusted for age, age2, sex, age*sex, age2*sex, 10 genetic principal components obtained from common genetic variants (MAF>1%), and 20 genetic principal components obtained from rare genetic variants (MAF<1%). The Regenie software [18] was used to perform all burden tests, and generate the scores above. Regenie uses Firth penalized likelihood to adjust for rare or unbalanced events, providing unbiased effect estimates.All analyses were performed separately for each of six genetic ancestries (African, Admixed American, East Asian, European, Middle Eastern, and South Asian). Summary statistics were meta-analyzed as for the single variant analysis. Participant assignment to genetic ancestry was done locally by each cohort, more details on the methods can be found in the S3 Table.Lastly, we used ACAT [35] to meta-analyze p-values across masks, within each phenotype separately. ACAT is not affected by lack of independence between tests. These values were used to draw Manhattan and QQ plots in Fig 2.Supporting informationS1 TableAcknowledgments per cohort.(XLSX)Click here for additional data file.S2 TableMembers of DeCOI, GEN-COVID (Italy), GEN-COVID (Spain), Mount Sinai Biobank, COVID-19 Host Genetics Initiative.(XLSX)Click here for additional data file.S3 TableCohorts demographics and sequencing details.(XLSX)Click here for additional data file.S4 TableInflation factors.(XLSX)Click here for additional data file.S5 TableSingle variant association tests significant results.(XLSX)Click here for additional data file.S6 TableFull Burden test results.(XLSX)Click here for additional data file.S7 TableIFN related genes burden test results.(XLSX)Click here for additional data file.S8 TableGWAS loci burden test results.(XLSX)Click here for additional data file.S9 TableAdditional information on the Genentech cohort.(XLSX)Click here for additional data file.S1 FigQQ plot and Manhattan plot for the exome-wide single variant association studies.(DOCX)Click here for additional data file.S2 FigChromosome 3 exome single variant association studies result by cohort for the severe disease phenotype.Note that all ancestries are shown together. Black dashed line represents the nominal statistical significance threshold (p = 0.05).(DOCX)Click here for additional data file.S3 FigChromosome 3 exome single variant association studies result by cohort for the severe hospitalized phenotype.Note that all ancestries are shown together. Black dashed line represents the nominal statistical significance threshold (p = 0.05)(DOCX)Click here for additional data file.S4 FigChromosome 3 exome single variant association studies result by cohort for the susceptibility phenotype.Note that all ancestries are shown together. Black dashed line represents the nominal statistical significance threshold (p = 0.05).(DOCX)Click here for additional data file.S5 FigSingle variant association study results at the three novel loci (OR and 95% CI).(DOCX)Click here for additional data file.S6 FigQQ plot from exome burden test ACAT meta-analyses.(DOCX)Click here for additional data file.S7 FigTLR7 beta coefficients (on logistic scale) with 95% interval for the severe COVID-19 phenotype, MAF<0.1%.X-axes are cut at -10 and 10.(DOCX)Click here for additional data file.S8 FigTLR7 beta coefficients (on logistic scale) with 95% interval for the severe COVID-19 phenotype, MAF<1%.X-axes are cut at -10 and 10.(DOCX)Click here for additional data file.S9 FigMARK1 beta coefficients (on logistic scale) with 95% interval for the severe COVID-19 phenotype, MAF<0.1%.X-axes are cut at -10 and 10.(DOCX)Click here for additional data file.S10 FigMARK1 beta coefficients (on logistic scale) with 95% interval for the severe COVID-19 phenotype, MAF<1%.X-axes are cut at -10 and 10.(DOCX)Click here for additional data file.S11 FigMARK1 beta coefficients (on logistic scale) with 95% interval for the hospitalized COVID-19 phenotype, MAF<0.1%.X-axes are cut at -10 and 10.(DOCX)Click here for additional data file.S12 FigMARK1 beta coefficients (on logistic scale) with 95% interval for the hospitalized COVID-19 phenotype, MAF<1%.X-axes are cut at -10 and 10.(DOCX)Click here for additional data file.S13 FigAncestry stratified results for TLR7, MARK1, and ABO.X-axis cut at 50. Figures show odds ratios and 95% confidence intervals.(DOCX)Click here for additional data file.
Medicine and Health Sciences Medical Conditions Infectious Diseases Viral Diseases Covid 19 Diagnostic Medicine Virus Testing Research and Analysis Methods Mathematical and Statistical Techniques Statistical Methods Metaanalysis Physical Sciences Mathematics Statistics Biology and Life Sciences Computational Biology Genome Analysis Genome-Wide Association Studies Genetics Genomics Human Genetics Biology and life sciences Organisms Viruses RNA viruses Coronaviruses SARS coronavirus SARS CoV 2 Microbiology Medical microbiology Microbial pathogens Viral pathogens Medicine and health sciences Pathology and laboratory medicine Pathogens Genetic Loci Genetics of Disease