Abstract
Importance Observational studies have suggested bidirectional associations between psychiatric disorders and COVID-19 phenotypes, but results of such studies are inconsistent. Mendelian Randomization (MR) may overcome limitations of observational studies, e.g. unmeasured confounding and uncertainties about cause and effect.
Objective To elucidate associations between neuropsychiatric disorders and COVID-19 susceptibility and severity.
Method In November, 2020, we applied a two-sample, bidirectional, univariable and multivariable MR design to genetic data from genome-wide association studies (GWASs) of neuropsychiatric disorders and COVID-19 phenotypes (released on 20 Oct. 2020). Our study population consisted of almost 2 million participants with either a (neuro)psychiatric disorder or data on COVID-19 status. Outcomes and exposures were anxiety, anxiety-and-stress related disorders, major depressive disorder, schizophrenia, bipolar disorder, schizophrenia-bipolar disorder combined (BIP-SCZ), and Alzheimer’s dementia on the one hand; and self-reported, confirmed, hospitalized, and very severe COVID-19 on the other.
Results In single-variable MR analysis the most significant and only Bonferroni-corrected significant result was found for BIP-SCZ (a combined anxiety of bipolar disorder and schizophrenia as cases vs. controls): the effect estimate was consistent with increased risk of COVID-19 (OR = 1.17, 95% CI, 1.06-1.28; p = 0.0012). Nominally significant univariable results were in line with slightly elevated risks of COVID-19 for genetic liabilities to bipolar disorder and schizophrenia. No COVID-19 phenotype consistently increased risk of (neuro)psychiatric disorders. In multivariable MR, bipolar disorder was the only phenotype showing a Bonferroni-corrected significant effect on a COVID-19 phenotype, namely severe COVID-19 (OR = 1.293; 95% CI, 1.095-1.527; p = 0.003). All sensitivity analyses confirmed the results.
Conclusions Genetic liability to bipolar disorder slightly increases COVID-19 susceptibility and severity. The contribution of bipolar disorder to these COVID-19 phenotypes was smaller than the odds ratios estimated by observational studies. Strength of association and direction of effect for genetic liability to schizophrenia were similar, albeit less significant. We found no consistent evidence of reverse effects, i.e. of genetic liability to COVID-19 on psychiatric disorders.
Introduction
Several large population-based have investigated associations between positive testing for COVID-19 on the one hand and psychiatric disorders on the other 1-3. Positive test result likelihoods for psychiatric disorders are inconsistent between those studies: while two of those cohort studies (from the UK and South Korea) do not report positive associations between COVID-19 testing and psychiatric disorders 1,4, others (from the UK and the US) mention odds ratios of 1.5 to 10 for associations between mental disorders and a COVID-19 diagnosis 2,3,5. For a recent diagnosis of a mental disorder in the US, odds ratios for COVID-19 were reported to be around 7.6, with evidence for relatively severe COVID-19 outcomes in those with a diagnosis of mental disorder5. One limitation of these studies is that psychiatric diagnoses were grouped by some 2,4, e.g. as mental disorders, psychotic disorders or affective disorders, precluding conclusions about COVID-19 risks for specific psychiatric disorders. Additionally, some diagnoses, such as bipolar disorder, were not included in some of the analyses 3. Furthermore, correction for medical comorbidities decreased several high odd ratios5. Finally, some authors admit residual socioeconomic factors may not be sufficiently captured in some databases2. Indeed, as recently noted, confounding or biases (partly) explain associations of COVID-19 with a range of traits and diseases 6. Mendelian Randomization (MR) has the potential to overcome two major limitations of observational studies: unmeasured confounding and uncertainties about cause and effect. Examples of MR studies elucidating risk factors for COVID-19 include two recent studies showing that BMI and smoking are risk factors for COVID-19 7,8. We are not aware of preprinted or published MR studies of psychiatric disorders and COVID-19. The recently updated whole-genome data on COVID-19 phenotypes (https://www.covid19hg.org/results/, see methods section) further increases the timeliness of MR approaches to elucidate risk factors for COVID-19 diagnosis and severity.
We hypothesized that given the aforementioned inconsistent observational evidence and attenuation of reported odds ratios when including covariates, psychiatric disorders do not constitute strong risk factors to contract COVID-19 or develop a severe course of COVID-19. Similarly, we hypothesized that COVID-19 would not increase susceptibility to psychiatric disorders. To test our hypotheses, we performed a range of bidirectional univariable and multivariable MR analyses of genetic liability to major psychiatric disorders and COVID-19 susceptibility as well as severe COVID-19.
Methods
Overview
We performed 2-sample mendelian randomization (MR) using summary statistics from large genome-wide association study (GWASs) to examine whether genetic liabilities to (neuro)psychiatric disorders increase risk of contracting COVID19 and of a severe course of COVID-19 (forward MR analyses, considering (neuro)psychiatric disorders as exposures and COVID-19 as outcome). In addition, we used COVID-19 GWASs to examine potential effects of genetic liabilities to COVID-19 diagnosis and to COVID-19 severity on (neuro)psychiatric diseases (reverse MR analyses, considering COVID-19 as the exposure and (neuro)psychiatric phenotypes as the outcomes). Most of our methods outlined below have been previously explained in more detail 9. The GWAS summary statistics we used were drawn from studies that had obtained written informed consent from participants and received ethical approval from institutional review boards. No ethical approval for the current analyses was needed as they were based on publicly available summary statistics.
(Neuro)psychiatric summary statistics
We used the available GWASs with summary statistics for psychiatric disorders including Alzheimer’s dementia (AD) (a meta-analysis of the stage 1 discovery dataset 10), anxiety 11,12, anxiety and stress-related diagnoses (ASRD) 11,12, Major Depressive Disorder (MDD) 13, Bipolar disorder (BIP) 14, Schizophrenia (SCZ) 15, and a meta-analysis of BIP and SCZ (BIP-SCZ; Table 1)16. To avoid that sample overlap between exposures’ and outcomes’ datasets impacted results substantially by inducing instrument bias in MR analyses, we excluded UK Biobank (UKBB) cohorts from (neuro)psychiatric disorders summary statistics. The largest Anorexia Nervosa GWAS and the largest 2018 and 2019 MDD GWASs contained a large number of UKBB participants, precluding us from using those GWASs for MR analyses.
GWASs of (neuro)psychiatric disorders used for the current study and outcomes used for the univariable forward MR analyses. In bold are depicted the (neuro)psychiatric GWASs that had identified ≥2 genome-wide significant loci and were thus selected as exposures in our forward MR analyses.
Hence, we selected MDD summary statistics from the 23andme cohort13, containing only 10,000 independent SNPs due to participants consent, and thus used these summary statistics only for univariable forward MR analysis. Similarly, we refrained from using the psychiatric cross-disorder GWAS as its study population was also partly composed of UKBB participants. For anxiety and ADRS, we performed meta-analysis in METAL 17 excluding UKBB participants: of anxiety using the iPSYCH (4,584 cases and 19,225 controls) and ANGST cohorts (7,016 cases and 14,745 controls) and of ASRD also using the iPSYCH cohort (4,584 cases and 19,225 controls) and ANGST cohorts (12,665 cases and 14,745 controls). We used the anxiety and ASRD summary statistics we had thus generated for our MR analyses.
COVID-19 summary statistics
We used the most recent COVID19 (GWAS) meta-analyses freeze V4 18 released on October 20, 2020 from the COVID-19 Host Genetics Initiative (https://www.covid19hg.org/results/) containing several COVID-19 phenotypes. This GWAS was based on a study population drawn from multiple cohorts: BioMe, FinnGen, Genes & Health, LifeLines Global Screening Array, LifeLines CytoSNP, Netherlands Twin Register, Partners Healthcare Biobank, and UKBB (Table 2). The majority of the included subjects were of European descent, with a small proportion of other ethnic backgrounds. Only variants with imputation quality > 0.6 had been retained and meta-analysis of individual studies had been performed with inverse variance weighting by the authors of the COVID-19 Host Genetics Initiative. To avoid that sample overlap between exposures’ and outcomes’ datasets impacted results substantially by inducing instrument bias in MR analyses, we only included COVID-19 summary statistics that had excluded the 23andme cohort. We divided the COVID-19 phenotypes (A-D, Table 2) into two categories, COVID-19 susceptibility and severity. We used two phenotypes to index COVID-19 susceptibility, namely C (defined by the COVID-19 Host Genetics Initiative as partial susceptibility) and D (self-reported COVID-19). Similarly, we used A (very severe respiratory confirmed COVID-19, here: ‘very severe COVID-19’) and B (hospitalized lab confirmed COVID-19, here ‘severe COVID-19’) to index COVID-19 severity.
Phenotype definition of the COVID-19 phenotypes by the authors of the COVID-19 Host Genetics Initiative used for the current study and outcomes used for the univariable reverse MR analyses. The numbers listed are for the GWASs conducted excluding the 23andme cohorts that were used for the current analyses. For further detail please see https://www.covid19hg.org/results/. In bold are listed the COVID-19 GWASs that had identified ≥2 genome-wide significant loci and were thus selected to be used as exposures in our reverse MR analyses.
MR analyses
Data were analyzed in November, 2020. As our main analysis, we performed univariable MR using GSMR (Generalized Summary Mendelian Randomization) 19 implemented in the Genome Complex Trait Analysis (GCTA) software 20. For our exposures, we first selected all relevant single-nucleotide variants (SNVs) identified in each GWASs as having reached a selection P-value threshold < 5 × 10−8 and being non-duplicate and uncorrelated (10 000 kilobase pairs apart and Linkage Disequilibrium (LD) R2 ≤.001). Instrument outliers were identified using HEIDI-outlier test (p <0.01), with the minimum number of instruments required for the GSMR analysis (nsnps_thresh) = 2. In harmonizing exposure and outcome data we removed palindromic SNVs with intermediate allele frequencies, and SNVs with minor-allele frequency (MAF) differences > 0.2 between exposure GWASs and outcome GWASs. We estimated the F parameter to evaluate instrument strength.
To interpret our results, we advise readers to take into account effect sizes and not focus on p-values. As a cut-off for statistical significance, we Bonferroni corrected two-sided p<0.05 for the number of tests performed in all analyses, i.e. univariable forward (n=35;), univariable reverse MR (n=16), and multivariable forward and reverse (see results section).
Then, Bonferroni-corrected significant results from GSMR were validated by applying as sensitivity analyses several alternative MR models, namely fixed-effect inverse variance– weighted (IVW), MR-Egger 21, weighted median-based regression methods 22, and MR-PRESSO 23. The harmonized input data from GSMR were used to perform by “TwoSampleMR” R packages 24. We used the MR-Egger intercept test, Cochran Q heterogeneity test, and MR pleiotropy residual sum and outlier (MR-PRESSO) test to evaluate potential inverse variance (IV) violations. We also performed leave-one-out analyses to examine whether any high-impact instruments possibly influenced MR results disproportionally.
At last, we conducted multivariable MR(MVMR) analyses 25 using the MendelianRandomization R package to examine which phenotypes remained risk factors taking into account pleiotropic effects among exposures. For MVMR analyses, we constructed instruments using SNVs in each of the GWASs meeting our single-variable MR selection criteria, as described above. We used the MVMR extension of the inverse-variance weighted MR method and MR-Egger method to correct for measured and unmeasured pleiotropy. In forward MVMR, we used as exposures SCZ, BIP and AD, while A1, A2, B1, B2, C1, C2, and D1 COVID-19 phenotypes (Table 1) were the outcomes. We excluded the following (neuro)psychiatric diseases as exposures in forward MVMR: BIP-SCZ 16 as it was highly correlated with SCZ and BIP (with most of the SNVs overlapping with either SCZ or BIP); MDD 13 as it did not have enough SNV information (summary statistics containing only 10,000 independent SNPs due to participants consent); anxiety and ASRD 11,12 as they did not have enough instrument variables, i.e. <2. In reverse MVMR analyses, we used the three COVID-19 A2, B2 and C2 phenotypes (Table 2) as exposures since they had the required number of ≥2 genome-wide significant loci, while outcomes were ASRD, AD, anxiety, BIP, SCZ, and BIP-SCZ. We removed duplicate and correlated SNVs (within 10 000 kilobase pairs; LD R2 ≥0.001), resulting in 8 SNPs as COVID-19 instruments and 119 SNPs as (neuro)psychiatric instruments. Statistical significance for MVMR was Bonferroni-corrected for the number of outcomes (see results section).
Data and code availability
We have made our code publicly available on Github (https://github.com/Bochao1/MR_PSY_COVID19). The datasets we accessed to perform our analyses may be found in the publications that listed as references for the datasets used.
Results
Overview
Five neuropsychiatric disorders (numbers of instruments: AD = 32, BIP = 6, MDD =5, SCZ = 137 and BIP-SCZ = 96; total study populations: AD = 63,926; BIP = 198,882; MDD = 307,354; SCZ = 105,318; and BIP-SCZ = 107,620) had enough (≥2) genome-wide loci to perform forward MR analyses (Table 1) on the seven COVID-19 phenotypes defined in Table 2. Conversely, three COVID-19 phenotypes (numbers of instruments: A2 = 8, B2 = 7 and C2 = 3; total study populations: A2 = 1,308,275, B2 = 969,689, C2 = 1,388,510; Table 2) had enough genome-wide loci to perform reverse MR analyses. As the instrument strength was strong (F-statistic in forward and reverse MR analyses ranging from 36.2 to 69.5), we did not find evidence of weak instrument bias.
Forward MR results of (neuro)psychiatric disorders as potential risk factors for COVID-19 phenotypes
We performed 35 univariable MR tests to examine the potential effects on COVID-19 of genetic liabilities to 5 psychiatric phenotypes with at least 2 significant loci identified. We thus Bonferroni corrected the significance threshold to 0.05/35 = 0.0014. Single-variable GSMR analysis showed that the top MR result was BIP-SCZ (combined GWAS of bipolar disorder and schizophrenia as cases vs. controls): the effect estimate was consistent with increased COVID-19 susceptibility (D1; N=96 instruments; OR, 1.165, 95% CI, 1.062-1.277; P = 0.0012; Figure 1). All four sensitivity MR analyses confirmed the direction of effect detected by GSMR (Table 3; Supplementary Table 1). IVW and MR-PRESSO also showed similar p-values to GSMR (Table 3). Although MR Egger was not significant (p= 0.247), overall horizontal pleiotropic effects were absent (for intercept of MR Egger: b= 0.015, p = 0.293). The I2GX statistics (I2GX =1.07%) of MR Egger were substantially smaller than 90%, suggesting substantial bias in the causal estimates due to uncertainty in the genetic associations, resulting in MR-Egger results not being reliable. The Cochran’s Q test in the fixed-effect IVW model (Q statistic = 9.395, p = 0.152) and MR Egger model (Q statistic = 7.278, p = 0.2) suggested that there was absence of heterogeneity in the instrumental variables, which may be the result of true causality rather than violation of instrument variable assumptions. Furthermore, the leave-one-out analysis (Figure 2A) showed that no SNVs altered the pooled IVW beta coefficient, confirming the stability of our results. We also found absence of directional pleiotropy and the instrument strength independent of direct effect (InSIDE) assumption to be satisfied (Figure 2B)21.
Forward MR results of BIP-SCZ as a risk factor for self-reported COVID-19 (D1), using 96 instrument variables.
Scatter plot of MR analyses using several models to examine causal relationships of BIP-SCZ genetic liability on self-reported COVID-19. The five models applied in the current manuscript are all depicted. Lines in black, green, brown, blue and purple represent results for fixed-effect IVW, weighted median, MR Egger, GSMR and MR-PRESSO models using 96 instruments. Neither GSMR nor MR-PRESSO identified any instrument outliers. Hence, the MR-PRESSO result was same as the IVW result, which was almost the same as the GSMR result, resulting in overlapping lines in the graph.
Leave-one-out analysis to evaluate whether any single instrumental variable was driving the causal association of BIP-SCZ with self-reported COVID-19 disproportionately. As can be appreciated from the graph, no genetic variant altered the pooled beta coefficient, indicating stability of our results.
Generalized funnel plot of univariable MR analysis of BIP-SCZ genetic liability effects on self-reported COVID-19 with first-order IVW and MR-Egger regression slopes to look for asymmetry as a sign of pleiotropy. This kind of graph plots the ratio estimate for each variant on the horizontal axis against its square-root precision (or weight) on the vertical axis. As can be appreciated from the plot, no evidence for asymmetry was detected, indicating absence of directional pleiotropy and the instrument strength independent of direct effect (InSIDE) assumption to be satisfied.
At a nominally significant p-value threshold (p<0.05), we detected causal effects of genetic liability to SCZ on very severe COVID (A2) and COVID-19 susceptibility (D1), of genetic liability to BIP on severe and very severe COVID-19 (A1 and B2), and of genetic liability to AD on COVID-19 susceptibility (C1, Supplementary table 2A). In line with the significant GSMR finding for BIP-SCZ, all of these nominally significant results had positive odds ratios for (neuro)psychiatric disorders on COVID-19 phenotypes (Supplementary table 2A).
Reverse MR results of COVID-19 phenotypes as potential risk factors of neuropsychiatric disorders
Three COVID-19 phenotypes met our predefined cut-off for inclusion as exposures in univariable MR analysis, namely severe and very severe COVID-19 (A2, B2) as well as COVID-19 susceptibility (C2, Table 2). MDD did not have enough overlapping SNPs to extract ≥2 instrument variables, resulting in 16 tests performed in reverse analyses (6 for A2, 6 for B2 and 4 for C2; Table 2; and Supplementary Table 2B). The p-value cut-off for significance was thus Bonferroni corrected to 0.05/16 < 0.0031. We found no results withstanding this multiple testing correction. Among most significant results were: genetic liability to very severe COVID-19 (A2) increasing risk of SCZ (GSMR OR, 1.036, 95% CI, 1.003-1.071; p = 0.031) and of BIP-SCZ (GSMR OR, 1.034, 95% CI, 1.002-1.068; p = 0.040); and genetic liability to severe COVID (B2) increasing risk of SCZ (GSMR OR, 1.046, 95% CI, 1.002-1.092; p = 0.041, Supplementary Table 2B).
Multivariable MR analyses
In forward MVMR, we examined the potential effects of genetic liabilities to three (neuro)psychiatric phenotypes jointly (AD, BIP, and SCZ; see methods) on seven COVID-19 phenotypes, resulting in a Bonferroni corrected (for the number of exposures) p-value threshold of p < 0.05 / 7 = 0.0071 (Supplementary Table 3A). Genetic liability to BIP showed a robust relationship with severe COVID-19 (B2; MVMR-IVW OR = 1.293; 95% CI, 1.095-1.527; p = 0.003). The estimates were consistent with estimates from the MVMR-Egger sensitivity analysis (OR = 1.302; 95% CI, 1.089-1.556; p = .005). The MR-Egger intercept analysis again did not indicate horizontal pleiotropy (Egger intercept: -0.015 CI, -0.046-0.016; p = 0.35; Supplementary Table 3A).
For reverse multivariable MR analyses examining the potential effects of genetic liabilities to three COVID-19 phenotypes jointly (A2, B2 and C2, i.e. GWASs with ≥2 genome-wide significant hits) on (neuro)psychiatric phenotypes, no COVID-19 phenotypes showed a causal relationship with any of neuropsychiatric disorders at a Bonferroni-corrected or nominal significance level in both IVW and MR-Egger (Supplementary Table 3B). The estimates were consistent with estimates from the MVMR-Egger sensitivity analyses (Supplementary Table 3B).
Discussion
We evaluated potential bidirectional associations between (neuro)psychiatric diseases and COVID-19 susceptibility and severity. Our univariable and multivariable MR results jointly indicate genetic liability to bipolar disorder slightly increases susceptibility to COVID-19 and COVID-19 severity. A combined phenotype of bipolar disorder and schizophrenia showed evidence in univariable MR for an increased risk of contracting COVID-19. This combined phenotype was not included in multivariable MR given substantial overlap with bipolar disorder and schizophrenia GWASs.
Our findings are consistent with observations in recent cohort studies in the USA and the UK, although the odds ratios we found were well below those reported in these studies 2,3,5. Contrary to one of these studies 2, we found no evidence of COVID-19 influencing risk of psychiatric disorders.
One possible explanation for the increased chances of contracting COVID-19 in those with a diagnosis of bipolar disorder may relate to disinhibition or impulsivity, personality traits that have been found to be associated with bipolar disorder26, even outside manic episodes27. Possibly, patients with bipolar disorder have on general more impulsive personality traits, may sometimes forget about social distancing and thus increase their chances of contracting COVID-19. Alternatively, manic episodes may play a more decisive role, increasing the odds of contracting COVID-19 for a variety of reasons (including risky behaviors) only during mania. Overall, however, the effect size was small, suggesting only slightly elevated risks of contracting COVID-19 and of COVID-19 severity when having genetic liability to bipolar disorder. Nonetheless, an implication of our finding may be awareness in patients as well as clinicians, resulting in attention to these elevated risks in outpatient consultations as well as emergency settings, e.g. those resulting from manic episodes.
Strengths of our study include the integration of univariable and multivariable, bidirectional MR analyses using many instrument variables drawn from large GWASs. We also included a range of (neuro)psychiatric as well as COVID-19 phenotypes and validated our results across available MR methods. Furthermore, in follow-up analyses such as leave-out analysis and using regression slopes of MR-Egger and IVW we demonstrate the robustness of our results. Finally, employing a multivariable in addition to a univariable approach is a strength, which was underscored by the converging findings for bipolar disorder. However, our results should also be interpreted in light of several limitations. First, a general concern in MR studies is risk of sample overlap. We minimized chances of sample overlap between exposure datasets and outcome datasets by excluding UKBB populations from (neuro)psychiatric GWASs and by excluding 23andme cohorts from COVID-19 datasets. Nonetheless, cryptic relatedness and potential sample overlap between exposure and outcome GWASs may result in some degree of instrument bias. However, the F-statistics we found were all above 36, allaying major concerns about weak instrument bias. Another limitation directly follows from the availability of GWAS data. For some phenotypes, such as obsessive-compulsive disorder and anorexia nervosa, no large datasets excluding the UKBB were available at the time of analysis or writing. For MDD and SCZ, summary statistics of larger GWASs may become available in the coming year. Similarly, as COVID-19 GWAS sample sizes ramp up, statistical power in MR analysis may increase. We encourage researchers to repeat MR analyses on other phenotypes and to use such larger GWAS datasets once they become available. To that end, we have uploaded our code to Github (see data availability section above).
In conclusion, we provide converging evidence for slightly increased susceptibility to and severity of COVID-19 in those with genetic liability to bipolar disorder. Odds ratios and direction of effect for genetic liability to schizophrenia were similar, albeit less significant. The contributions of bipolar disorder and schizophrenia to these COVID-19 phenotypes were smaller than the odds ratios estimated by observational studies. We found no consistent evidence of reverse effects, i.e. of genetic liability to COVID-19 on psychiatric disorders.
Data Availability
Data and code availability We have made our code publicly available on Github (https://github.com/Bochao1/MR_PSY_COVID19). The datasets we accessed to perform our analyses may be found in the publications that listed as references for the datasets used.
Conflict of interest statement
The authors declare no conflict of interest.
Funding
No funding was provided to carry out this project.