Genetic and phenotypic analysis of the causal relationship between aging and COVID-19.

Background: Epidemiological studies revealed that the elderly and those with comorbidities are most affected by COVID-19, but it is important to investigate shared genetic mechanisms between COVID-19 risk and aging.
Methods: We conducted a multi-instrument Mendelian Randomization analysis of multiple lifespan-related traits and COVID-19. Aging clock models were applied to the subjects with different COVID-19 conditions in the UK-Biobank cohort. We performed a bivariate genomic scan for age-related COVID-19 and Mendelian Randomization analysis of 389 immune cell traits to investigate their effect on lifespan and COVID-19 risk.
Results: We show that the genetic variation that supports longer life is significantly associated with the lower risk of COVID-19 infection and hospitalization. The odds ratio is 0.31 (P = 9.7 × 10-6) and 0.46 (P = 3.3 × 10-4), respectively, per additional 10 years of life. We detect an association between biological age acceleration and future incidence and severity of COVID-19 infection. Genetic profiling of age-related COVID-19 infection indicates key contributions of Notch signaling and immune system development. We reveal a negative correlation between the effects of immune cell traits on lifespan and COVID-19 risk. We find that lower B-cell CD19 levels are indicative of an increased risk of COVID-19 and decreased life expectancy, which is further validated by COVID-19 clinical data. Conclusions: Our analysis suggests that the factors that accelerate aging lead to an increased COVID-19 risk and point to the importance of Notch signaling and B cells in both. Interventions that target these factors to reduce biological age may reduce the risk of COVID-19.

Introduction

The coronavirus disease 2019 (COVID-19), caused by severe acute respiratory coronavirus 2 (SARS-CoV-2), first emerged in late 2019 and has led to an unprecedented global health crisis[1]. Notably, the aging population is particularly at risk of COVID-19[2], e.g., in Italy, 88% of the individuals who tested positive for COVID-19 were 40 years or older[3]. A recent report based on epidemiological data from multiple countries showed that 69% of infections in people over 70 are symptomatic, whereas this number drops to 21% for 10−19-year-olds[4]. Unsurprisingly, elderly people are also more likely to die from COVID-19, and the case fatality rate for COVID-19 grows exponentially with age[3]. As observational evidence implies a strong link between COVID-19 and age, COVID-19 can be considered a disease of aging[3], and multiple clinical trials using potential lifespan-extending drugs (e.g., metformin, rapamycin, and senolytics) to protect the elderly from COVID-19 have been proposed[5-7]. Although observational data on metformin seems promising[8,9], it is unclear if other lifespan-extending drugs should be prioritized in clinical trials since the evidence of any causal link between lifespan and COVID-19 susceptibility is still missing.

Mendelian Randomization (MR) is a genetic instrumental variable approach that assesses the causal effect of exposure of interest on an outcome by ascertaining genetic variants, e.g., single nucleotide polymorphisms (SNPs), strongly associated with the exposure phenotype. Since the alleles of the genetic variants are naturally randomly allocated at conception, when the genetic effects on the outcome are only mediated through the exposure, the causal effect inferred by MR is, in analogy to randomized clinical trials (RCTs), free of any environmental confounders and reverse causation. Although RCTs are considered a gold standard for establishing causal relationships, MR can provide valuable insights into causality when it is not feasible to perform an RCT or before an RCT is performed[10].

In this study, we perform a multi-SNP MR analysis to elucidate whether and how the rate of aging is associated with COVID-19. We consider four lifespan-related traits (parental lifespan, healthspan, longevity, and healthy aging (the combination of these three traits)) as exposures and evaluate their causal effects on COVID-19 infection and related phenotypes. To support the argument, we also estimate the biological age acceleration in COVID-19 patients from the UK Biobank (UKBB) cohort and observe a significant association between the phenotypic indicators of aging progress (aging clocks) and the risk and case fatality rate of COVID-19. To provide functional insight into how aging contributes to a higher risk of COVID-19, we conduct a bivariate genomic scan to highlight the loci contributing to both aging and COVID-19 risk, identifying the Notch signaling pathway and immune system development. Finally, we perform MR using 389 immune cell traits as exposure and observe a significant negative correlation between their effect on lifespan and COVID-19 risk, especially for B cell-related traits. More specifically, we discover that lower CD19 levels in B cells may increase the risk of COVID-19 and decrease lifespan, which is further validated by clinical data from COVID-19 subjects.

Methods

GWAS data for lifespan-related traits and diseases

We studied four lifespan-related traits (lifespan, longevity, healthspan, and a combined trait) with publicly available GWAS summary statistics. The parental lifespan GWAS included unrelated, European-ancestry subjects (a total of 512,047 mothers’ and 500,193 fathers’ lifespan), 60% of which were complete. The statistics for every cohort were calculated by fitting Cox survival models to mother’s and father’s survival, respectively, taking account of 10 principal components, study-specific covariates, and sex. In the GWAS setting, parental lifespan is the same phenotype as the general lifespan of individuals (but with a weaker power) due to the fact that the genetic effect on a parental phenotype is simply half of the individual’s phenotype itself. Thanks to the large sample size of UK Biobank, such a GWAS is powerful enough to uncover the genetic architecture[11].

The longevity GWAS included unrelated, European-ancestry subjects with a lifespan above the 90th survival percentile (N = 11,262) or whose age at the last follow-up visit (or age at death) was before the 60th percentile age (N = 25,483). The statistics for each cohort were calculated using logistic regression and then combined using a fixed-effect meta-analysis[12]. The healthspan GWAS contained 300,477 unrelated, British-ancestry individuals from UKBB. The statistics were calculated by fitting Cox−Gompertz survival models. The healthspan was defined as the age of the first incidence of dementia, congestive heart failure, diabetes, chronic obstructive pulmonary disease, stroke, cancer, myocardial infarction, or demise[13].

The summary association statistics of healthy aging was from the meta-analysis of healthspan, lifespan, and longevity summary statistics using MANOVA[14], while accounting for correlations between studies due to sample overlap and correlation amongst the traits. Summary association statistics were calculated for 7,320,282 SNPs shared between the studies. These statistics represented the significance of each SNP affecting one or more of the traits, giving a P-value against the null hypothesis that effect sizes are zero in all studies[14,15].

We investigated four additional traits genetically correlated with lifespan using published case-control studies: Alzheimer’s disease[16], coronary artery disease (CAD)[17], type 2 diabetes[18], and smoking[19] (Table S1).

We also included GWAS for age acceleration measured by four epigenetic clocks, including Hannum age, Horvath age, PhenoAge, and GrimAge[20]. The epigenetic age was estimated for 34,449 healthy individuals of European ancestry. In addition to epigenetic age, we include two physical function-related traits, the pace of walk and the sedentary lifestyle, as they are correlated with the rate of aging and therefore can serve as the surrogates to the biological age[21,22].

GWAS data for 22 common diseases were from a community-based study, Genetic Epidemiology Research on Adult Health and Aging (GERA)[23]. There were 60,586 individuals of European ancestry in the GERA data. The summary statistics of these diseases were adjusted with age, gender, and the first 20 PCs.

We used 1000 Genomes Phase 3 reference (released in 2014 October) to map variants in the GWAS results to rsIDs by chromosome, position, and alleles. Only the autosomal SNPs available in the 1000 Genomes reference panel were used, and the 1000 Genomes European ancestry reference was used to estimate the linkage disequilibrium (LD) among these SNPs. Duplicated rsIDs in the data were removed prior to the analysis.

COVID-19-related traits

To extensively evaluate the genetic effects on COVID-19 risk, we used GWAS summary statistics data from 8 COVID-19-related traits (Table S1). The GWAS results for SARS-COV-2 infection are from the National Institute of Health, Genome-Wide Repository of Associations Between SNPs and Phenotypes (NIH-GRASP), which includes 1,503 positive cases and 11,409 negative or 457,747 UK Biobank controls with European ancestry; the GWAS summary statistics for the critical illness was from the GenOMICC (Genetics Of Mortality In Critical Care) study in 2,244 critically ill Covid-19 patients from 208 UK intensive care units[24]. The rest of the five traits are from the COVID-19 Host Genetics Initiative (HGI) release 5 (Jan 2021), with the sample size varies from 1,332 to 1,079,768[25]. Those traits including COVID-19 hospitalization (versus non-hospitalized COVID-19 or population control), susceptibility (affected versus unaffected population), very severe respiratory confirmed COVID-19 (versus the general population), and COVID-19 infection (versus population).

Expression quantitative trait loci (eQTLs) and age-related gene expression in blood

Blood eQTL data were obtained from the eQTLGen Consortium (31,684 whole blood samples)[26]. Only the significant near-independent eQTLs (FDR-q < 0.05, r2 < 0.05) were used in the MR analysis.

The age-related transcriptomic change in whole blood was obtained from a large-scale meta-analysis[27], including six European-ancestry studies (n = 7,074 samples), and detected roughly half of the genes in the human genome (n = 11,908). The direction and P-value of age-related differential expression were directly obtained from the published dataset.

Immune cell traits

The GWAS summary statistics of immune cell-type-specific surface marker levels are obtained from the largest immune cell GWAS study[28]. 389 median fluorescence intensities (MFIs) of surface antigens were profiled by flow cytometry and assessed in a general population cohort of 3,757 Sardinians.

Mendelian randomization analysis

MR is a method that uses genetic variants as instrumental variables to determine whether an observational association between a risk factor and an outcome is consistent with a potential causal effect[29]. The multi-SNP MR analysis was implemented using GSMR (Generalized Summary-data-based MR) in GCTA[30].

As instruments for each exposure (four lifespan-related traits, four risk factors, and four epigenetic age acceleration traits), we selected near-independent SNPs (r2 < 0.1) with genome-wide significant (P < 5 × 10−8) association with the exposure. For the expression of NOTCH1-4 in whole blood and other tissues, we selected significant near-independent eQTLs (FDR-q < 0.05, r2 < 0.05); For 22 diseases from GERA community-based study, we selected SNPs with suggestive genome-wide significance (P < 1 × 10−6) as instruments and performed a separate analysis due to the limited case number in the community-based study. A full list of genetic instruments is provided (Supplementary Data 1).

GSMR includes a HEIDI-outlier filter to remove potential pleiotropic SNPs that affect the exposures and the outcomes via different pathways. We set its p-value threshold to 0.01 and tested the remaining SNPs for association with the COVID-19-related traits. The required minimum number of instrumental SNPs for each exposure in the analysis was lowered to 1.

Conditional analysis

To test whether the effect of lifespan-related traits on COVID-19 risk depends on certain age-related diseases and vice versa, we performed a conditional analysis using a two-step approach, as described by Zhu et al.[30]. In the first step, we performed a conditional GWAS analysis to adjust the exposure of interest with other risk factors using mtCOJO (multi-trait-based conditional and joint analysis). In the second step, we estimate the effect of adjusted exposure on the outcome using GSMR as previously described. We, therefore, can estimate the effect size of lifespan-related traits on COVID-19, accounting for other age-related risk factors by a GSMR analysis using SNP effects conditioning on covariate traits. Notably, as the exposures are very highly correlated, the multivariate MR will have lower power. To estimate the causal effects of conditional traits, we had to lower the P-value threshold for genetic instruments to 5e−6 to obtain a sufficient number of SNPs for MR analysis. To make the univariate and conditional analysis results comparable, we also performed a univariate analysis using the same P-value threshold.

Sensitivity analysis

We used GSMR for the main analyses because it gains power by taking account of sampling variation of the effect size of SNPs on exposure and outcome, compared with the MR-Egger and inverse variance weighted (IVW) methods[30]. GSMR also accounts for the remaining LD among instruments after clumping analyses. To compare the results from other MR methods based on various assumptions, we performed a sensitivity analysis using the Maximum likelihood method[31], MR-Egger method[32], and simple median method[33].

The Maximum likelihood method estimates the causal effect by maximization of the likelihood based on the effect of SNPs on exposure and outcome. It gives robust estimates even in the presence of small measurement error for the effect of SNPs on exposure[31]; the MR-Egger method is the modification of the IVW method, which allows a non-zero intercept. This way, it allows unbalanced pleiotropic effects across all of the instruments while still returns unbiased causal effect estimates. This method assumes no correlations between horizontal pleiotropic effects and SNP-exposure effects (the InSIDE assumption)[32]. The MR-Egger regression also provides an intercept test to detect the directional pleiotropy in the instruments (i.e., the pleiotropic effect is evident if the intercept term significantly deviates from 0). Lastly, the simple median method takes the median effect of all instrumental SNPs. It only requires half of the SNPs to be valid to return accurate causal effect estimates.

We selected independent instrumental SNPs (r2 < 0.01) for each exposure with the same genome-wide significance threshold as in GSMR analysis (P < 5 × 10−8 for lifespan-related traits and FDR < 0.05 for eQTLs); the analysis was then performed using the “TwoSampleMR” R package (https://mrcieu.github.io/TwoSampleMR)[34].

Bivariate genomic scan and functional annotation

To identify genetic variants associated with aging-related COVID-19 risk, we meta-analyzed UKBB COVID-19 infection (with population control) and healthy aging (with the sign of effect size reversed) summary statistics while accounting for correlations between studies due to sample overlap and correlation between the traits, as implemented in MultiABEL v1.1-610[14,35]. Summary association statistics were calculated for the 7,318,649 SNPs shared between studies. These statistics represent the significance and consistency of each SNP affecting one or both of the traits (e.g., the SNPs that significantly contribute to aging and COVID-19 risk in the same direction will have a smaller P-value). Therefore, we refer to this bivariate genomic scan result as the aging-related COVID-19 throughout this study.

We then used the summary statistics of aging-related COVID-19 and performed functional annotation for all SNPs in genomic areas identified by lead SNPs (P < 1 × 10−6, 250 Kb apart) using FUMA (Functional Mapping and Annotation)[36]. The annotated genes were used for functional enrichment analyses using the default setting of the FUMA platform.

Genetic correlation analysis

We estimated genetic correlations between lifespan-related traits, risk factors, epigenetic age acceleration, and COVID-19 using LD score regression (LDSC) and high-definition likelihood (HDL) methods[37,38]. SNPs that are imperfectly imputed (INFO < 0.9) or with low frequency (MAF < 0.05) were removed to reduce statistical noise. LDSC was performed using LDSC software v1.0.1 (https://github.com/bulik/ldsc); the HDL was performed using R package “HDL” v1.3.8 (https://github.com/zhenin/HDL).

Biological age estimation for UKBB cohorts

The collection of the UK Biobank (UKBB) data was approved by the UKBB’s Research Ethics Committee. Access to the UK Biobank data was granted for this work under UK Biobank application number 21988. All-cause mortality increases exponentially with age. Hence, log-linear risk predictors from proportional hazards models can provide natural composite organism state representations characterizing the progression of aging based on biological and physiological measurements. We used two such biological age measures: Phenotypic Age based on blood biochemistry[39] and Dynamic Organism State Indicator (DOSI) based on widely available Complete Blood Counts (CBC)[40]. The latter is a proxy for the frailty index and is derived from the blood markers only, whereas the Phenotypic Age additionally employs the explicit age. We also used physical activity (number of steps per day averaged over the week), which is associated with all-cause mortality and hence can also be viewed as a measure of biological aging[41].

We investigated an association between the incidence of COVID-19 and biological age acceleration (which is the difference between the biological age of a person and the average biological age in the cohort of individuals of the same age and sex) using logistic regression. Chronological age and biological sex were used as additional covariates in the analysis.

Following UKBB recommendations, we used the “result” label from the table “COVID-19 test results table” as the proxy of disease severity. This implies that primarily those individuals that showed characteristic COVID-19 symptoms were selected for testing. We investigated biological age acceleration associations with the incidence of COVID-19 and its associated fatality using all available cases (All) and separately cohorts of individuals who have (Frail) or do not have (Not Frail) major chronic diseases (from the list including all kinds of cancer, angina pectoris, coronary heart disease, heart attack, heart failure, hypertension, stroke, diabetes, arthritis, bronchitis, and emphysema) at the time of infection.

Table 1

Mendelian randomization estimates for the association between lifespan-related traits and risk of COVID-19.

Exposure	Outcome	OR	95% CI	P
Healthy aging	HGI covid susceptibility	0.33	0.13−0.85	2.2e−02
Healthy aging	HGI hosp covid vs. nonhosp	0.38	0.18−0.78	9.1e−03
Healthy aging	UKBB covid vs. neg	0.25	0.09−0.69	7.0e−03
Healthy aging	UKBB covid vs. pop	0.12	0.05−0.32	1.6e−05
Longevity	HGI covid susceptibility	0.68	0.56−0.82	8.5e−05
Longevity	HGI hosp covid vs. nonhosp	0.81	0.71−0.93	2.3e−03
Longevity	HGI hosp covid vs. pop	0.89	0.82−0.96	2.9e−03
Longevity	UKBB covid vs. neg	0.58	0.47−0.72	5.1e−07
Longevity	UKBB covid vs. pop	0.53	0.43−0.65	2.3e−09
Lifespan	HGI covid susceptibility	0.45	0.27−0.77	3.2e−03
Lifespan	HGI hosp covid vs. nonhosp	0.46	0.3−0.71	3.3e−04
Lifespan	HGI hosp covid vs. pop	0.71	0.55−0.91	6.8e−03
Lifespan	UKBB covid vs. neg	0.44	0.26−0.77	3.6e−03
Lifespan	UKBB covid vs. pop	0.31	0.18−0.52	9.7e−06

Only the associations that reached nominal significance (P < 0.05) are shown. hosp: hospitalized COVID-19 patient; nonhosp: non-hospitalized COVID-19 patient; pop: population control; neg: COVID-19 negative control.

Table 2

Association between biological age acceleration and the risk of COVID-19.

Biological age measurement	Outcome	OR	95% CI	P
Phenotypic age	COVID19 incidence (All)	1.28	1.25−1.31	8.4e−82
Phenotypic age	Case fatality (All)	1.19	1.04−1.35	1.1e−02
Phenotypic age	COVID19 incidence (Not frail)	1.12	1.04−1.2	1.9e−03
Phenotypic age	Case fatality (Not frail)	1.72	1.17−2.51	5.4e−03
Phenotypic age	COVID19 incidence (Frail)	1.26	1.23−1.3	3.7e−62
DOSI	COVID19 incidence (All)	1.31	1.26−1.38	9.5e−32
DOSI	COVID19 incidence (Not frail)	1.09	1.01−1.19	3.6e−02
DOSI	Case fatality (Not frail)	2.44	1.45−4.06	7.7e−04
DOSI	COVID19 incidence (Frail)	1.35	1.28−1.42	3.6e−28
Physical activity	COVID19 incidence (All)	0.95	0.93−0.96	9.1e−19

Only the associations that reached nominal significance (P < 0.05) are shown.

Table 3

Mendelian randomization estimates for the association between NOTCH2 expression and risk of COVID-19.

Exposure	Outcome	OR	95% CI	P
NOTCH1	Covid critical illness	2.57	1.39−4.74	0.0025
NOTCH2	UKBB covid vs. neg	1.46	1.08−1.99	0.0150
NOTCH2	UKBB covid vs. pop	1.43	1.07−1.91	0.0150

Only the associations that reached nominal significance (P < 0.05) are shown.