Introduction

The infectious disease syndrome known as coronavirus disease 2019 (COVID-19) is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and has become pandemic in 2020. Although most COVID-19 cases result in mild symptoms, some patients suffer severe symptoms, including severe pneumonia and multiorgan failure. There is also great variability in the duration of COVID-19 symptoms, with a majority of patients reporting remission of symptoms within 14 to 21 days in the outpatient setting1, while some patients experience prolonged multiorgan dysfunction and complications that last for 12 weeks or longer2. The study of host genetics can bring valuable support toward a better understanding of the mechanisms underlying COVID-19 and can guide the development of preventive and therapeutic measures to mitigate the health burden of this pandemic. Multiple efforts are underway to study the contribution of host genetics to COVID-19, with a focus on risk of severe COVID-19 outcomes and risk of infection with the SARS-CoV-2 virus3,4,5,6,7. In particular, one genetic locus, 3p21.31, has repeatedly been associated with severe respiratory illness and hospitalisation for reasons of COVID-195,7. However, few studies have focused on the genetics of symptoms duration and remission in outpatients. Persistence of COVID-19 symptoms is associated with a higher risk of complications, including prolonged hospitalization, and poor quality of life thereafter8.

Recently, the COLCORONA randomised clinical trial compared the benefit of low-dose colchicine to placebo in 4488 outpatient individuals diagnosed with a COVID-19 infection who were 40 years or older and with at least one high-risk criterion for severe disease9. The composite primary endpoint of death or hospitalisation for reasons of COVID-19 occurred in 4.7% of patients in the colchicine group and 5.8% of those in the placebo group (odds ratio (OR), 0.79; 95% confidence interval (CI) 0.61–1.03; P = 0.08). In a prespecified analysis of 4159 participants who received a diagnosis of COVID-19 confirmed by a polymerase chain reaction (PCR) test, the primary endpoint occurred in 4.6% and 6.0% of patients in the colchicine and placebo groups respectively (OR 0.75; 95% CI 0.57 to 0.99; P = 0.04)9. Participants were invited to take part in an optional genetic substudy to this randomized clinical trial. Here, we conducted a genetic study of time to remission of COVID-19 symptoms assessed during the 30-day follow up period of the COLCORONA clinical trial, with the aim to gain a better understanding of the underlying factors responsible for the duration of symptoms in the acute phase of disease in a recently diagnosed outpatient population.

Methods

Study population

COLCORONA was a randomised, double-blind, placebo-controlled trial comparing colchicine (0.5 mg twice daily for the first 3 days and then once daily for 27 days thereafter) with placebo in a 1:1 ratio10. Patients were eligible if they were at least 40 years of age, had received a diagnosis of COVID-19, were not hospitalised, and presented at least one of the following high-risk criteria: age of 70 years or more, obesity (body-mass index of 30 kg/m2 or more), diabetes, uncontrolled hypertension with systolic blood pressure > 150 mmHg, known respiratory disease, known heart failure, known coronary disease, fever of at least 38.4 °C within the last 48 h, dyspnea at the time of presentation, bicytopenia, pancytopenia, or the combination of high neutrophil and low lymphocyte counts. Exclusion criteria included inflammatory bowel disease, chronic diarrhea or malabsorption, pre-existent progressive neuromuscular disease, estimated glomerular filtration rate (eGFR) less than 30 mL/m in/1.73 m2, cirrhosis, chronic active hepatitis or severe hepatic disease, female patients who are pregnant or breast-feeding or are considering becoming pregnant during the study or for 6 months after the last dose of study medication, history of an allergic reaction or significant sensitivity to colchicine, concomitant chemotherapy for cancer, and patients who were considered by the investigators of COLCORONA trial, for any reason, to be an unsuitable candidate for the study were excluded from the original study. Clinical evaluation visits occurred by telephone at 15- and 30-days following randomisation from March 2020 to January 2021.

At the time of consent to the main COLCORONA study, 3315 participants were asked whether they could be recontacted to take part in the optional pharmacogenomic substudy. Those who agreed were called back by the team to obtain a separate consent for the genetic study and to plan for saliva sample collection. Exclusion details are shown in Supplementary Fig. 1. There were 2441 participants recruited into the genetic study. Of the 2249 samples genotyped, 20 were excluded due to < 98% genotyping completion rate, 3 samples with discordant sex between clinical and genetic data sets were excluded, one contaminated sample was removed, 16 genetically-determined related family members were excluded and we excluded 232 outliers from the cluster according to multidimensional scaling based on 1000 Genomes CEU reference samples identified as “Utah Residents with Northern and Western European Ancestry”. One patient was excluded from the intent to treat population, 3 never received the study drug, and 142 patients who did not have a confirmatory COVID-19 diagnostic test were excluded, leaving 1855 patients for analysis. Written informed consent was obtained from all participants. The study was approved by the Montreal Heart Institute research ethics committee and complies with the Declaration of Helsinki.

Endpoint definitions

The day of the end of COVID-19 symptoms was reported by the COLCORONA trial participants during the follow-up telephone visits that occurred on days 15 and 30 of the study. The event date was set as the study day of last known symptoms, alternatively patients were censored at the date of the patients’ last telephone visit. Hospitalized participants were considered as having ongoing symptoms. There were 132 (7%) participants with missing or incomplete information on reported symptom dates, leaving 1723 participants for the analysis of time to symptom remission. The proportion of participants with missing information in the genetic substudy was comparable to the proportion in the COLCORONA trial (359/4488; 8%). Further investigations for the 132 excluded participants found related research nurses’ notes for 34 patients, 14 had reported symptom initiation after the baseline visit, and 20 had reported having no symptoms. Information on the dates of symptom occurrence was collected, but the specific symptoms themselves were not collected. No deaths were reported in participants to the genetic study, consequently the severity outcome was limited to hospitalisation for reasons of COVID-19 in the 30 days following randomisation.

Genotyping and imputation

DNA was obtained using the Oragene OG-500 self-collection saliva kits delivered to and collected from the participants’ home while respecting distancing and quarantine restrictions. DNA was extracted from 1 ml of saliva and genome-wide genotyping was performed using 200 ng of genomic DNA at the Beaulieu-Saucier Pharmacogenomics Centre (Montreal, Canada). The Illumina Infinium Global Screening Array (GSA) v3-MD (Illumina, San Diego, CA) including 700,625 genomic markers was used and processed according to the manufacturer's specifications. BeadChips were scanned using the Illumina iScan Reader and analysed with the data manifest MHI_GSAMD-24v3-0-EA_20034606_C2.bpm, with minor manual cluster adjustment for ADME genes. Plink files were produced by the iaap-cli tool (version 1.1.0-80d7e5b). Intensities, B allele frequency, and log R ratio were extracted using the gtc_convert tool (version 0.1.2). pyGenClean11 version 1.8.3 and PLINK12 version 1.07 and 1.9 (the latter for the data manipulation steps of the relatedness and ethnicity modules) were used to conduct quality control and genetic data clean-up processes. The genotyping experiment consisted of 28 plates of DNA samples with one control per hybridization experiment selected from NA18861, NA06994, NA19147 and NA12878 obtained from the NIGMS Human Genetic Cell Repository at the Coriell Institute for Medical Research. The pairwise concordance of control samples ranged from 0.99998 to 0.999998 and concordance of their genotypes with expectation from the 1000 Genomes genotype data ranged from 0.9986 to 0.9989.

The completion rate threshold for genotypes and samples was set to 98%. SNPs with genotyping plate bias (n = 213) based on the 96 well plates used to dilute DNA samples were flagged but not removed as the effect of genetic ancestry could not be excluded. Genetic variants not in Hardy–Weinberg equilibrium in the European ancestry subgroup were excluded (P < 8.99 × 10–8). Pairwise identity-by-state (IBS) was used to conduct close familial relationship checks. We removed all but one member of related samples based on a selection of uncorrelated SNPs (r2 < 0.1). The pairwise IBS matrix was used as a distance metric to identify cryptic relatedness among samples and to identify cluster outliers by multidimensional scaling (MDS). The coordinates of the first two MDS components of each subject were plotted including the genotypes of HapMap CEU, JPT-CHB, and YRI data (unrelated individuals). Outliers from the main cluster overlapping the CEU reference samples (Utah residents with Northern and Western European ancestry from the CEPH collection) were removed according to k-nearest neighbour with a threshold of 1.9σ in pyGenClean (v1.8.3). Principal components were generated based on the final sample selection and used to control for confounding by ancestry in genetic association analyses13.

Genome-wide imputation was performed on the TOPMed Imputation Server (version 1.5.7)14 using Eagle (version 2.4)15 for phasing and Minimac4 (version 1.0.2)14 for imputation. Following quality control, 513,036 genetic variants were used for imputation. The imputation provided 308,070,060 (513,036 genotyped and 307,557,024 imputed) genetic variants from which we retained 43,062,846 (513,012 genotyped and 42,549,834 imputed) variants with a quality value (r2) ≥ 0.6, and of which 6,490,603 had a minor allele frequency (MAF) ≥ 5%. Genotype “hard calls” (non-probabilistic best-guess genotype assignations) were generated for individual genotypes with genotype probability scores ≥ 0.80 and otherwise set to missing. The hard calls were used only in analyses stratified by genotype. Chromosomal positions are according to GRCh38.

Statistical analyses

Genome-wide association analysis (GWAS) with common genetic variants (MAF ≥ 5%) were conducted using the program genetest version 0.5.016. Cox proportional hazards regression was used for the GWAS of time from randomisation to end of COVID-19 symptoms (time to remission) in 1723 study participants. Logistic regression was used for the analysis of the severity endpoint of hospitalisation for COVID-19 in 1855 participants. To avoid bias due to population structure, analyses were limited to unrelated individuals of genetically determined European ancestry, which was the largest population group in the sample. Genetic variants were coded according to imputation dosage in (0–2) and tested using a 1-degree of freedom Wald test adjusted for age, sex, and the first 10 principal components to control for genetic ancestry, with the addition of the study treatment arm when both arms were included in the analyses. GWAS were conducted using participants from both study arms, as well as stratified by treatment arms. We further assessed the statistical interaction between genetic variants and treatment arms to highlight possible treatment-specific effects. Sensitivity analyses were conducted with the two genetic variants identified by GWAS by adjusting for BMI, height, obesity, diabetes, hypertension, smoking status, history of respiratory disease, and coronary artery disease in addition to adjustment for age, sex, and the 10 principal components of genetic ancestry. The two genetic variants were also tested for association with height using a linear regression with adjustment for age and sex. Exploratory sex-stratified analyses were performed. The GWAS for hospitalisation due to COVID-19 was only performed using both study arms due to the limited number of events. Each GWAS was conducted at the 5 × 10–8 significance level to adjust for the multiple testing of genetic variants within each phenotype. No additional adjustments were made to account for the multiple phenotypes or different statistical methods. Results are reported with point estimates and 95% confidence intervals (CI) which are not adjusted for multiple comparisons. The top findings from the GWAS were reproduced using SAS version 9.4. The proportionality of hazards assumption was verified for variants identified in a Cox proportional hazards regression GWAS, and the convergence criterion was met for each considered model. Cumulative incidence plots were obtained using a parametric model fitted with the R package. In order to quantify the effect of the genetic variants upon the duration of symptoms, we fitted a Weibull accelerated failure time model adjusted for age, sex and 10 principal components using the lifereg procedure in SAS 9.4. To assess the representativity of the proportion of events in the genetic substudy to that of the main COLCORONA trial, we used the Fisher Exact test procedure to compare the number of patients with and without events in the genetic substudy to patients not in the genetic substudy.

Functional annotation

We first defined credible candidate variants as those located within 500 kb of the leading variants and with P values within two orders of magnitude of the lead variant. We used the software GCTA-CoLo17 to conduct a conditional analysis to identify independent signals. We used PAINTOR18 to identify credible sets of causal variants based on the magnitude and direction of association and the pairwise linkage disequilibrium structure at the loci, and we used RegulomeDB19 and DSNetwork20 to assign a relative ranking to variants. We used in silico functional annotations from the public databases Open Target Genetics21 and PhenoScanner22,23 to identify potential functional mechanisms and target genes. We tested the colocalization between the COLCORONA GWAS signals and clinically relevant phenotypes using the COLOC R package v3.2-124. Network analysis of selected candidate genes was conducted using GeneMANIA25. The complete approach is detailed in Supplementary Material.

Data availability

The anonymized patient level data from the COLCORONA trial will be shared via the Vivli (vivli.org) data repository. The patient level genetic data underlying this article cannot be shared to preserve the privacy of study participants. Summary statistics from the GWAS results are available for download and visualisation via PheWeb26 at statgen.org/pheweb/colcorona.

Results

There were 1855 participants available for the genetic study of COLCORONA (Supplementary Fig. 1). The baseline characteristics of the participants to the genetic study are shown in Table 1. The mean age of participants was 54.1 years, 56.2% were female, mean BMI was 30.2 kg/m2, 16.3% had a medical history of diabetes, 32.1% of hypertension, and 29.4% of respiratory disease. There were 1254 (72.8%) participants who reported being free of COVID-19 symptoms during the study follow-up period, and the mean number of days between randomisation and end of COVID-19 symptoms in those participants was 11.7 days. Mean number of days between first symptoms and randomisation was 5.4 days (Table 1). Overall, women reported longer symptom duration than men (Supplementary Table 1). The COLCORONA primary endpoint of death or hospitalisation for COVID-19 differed significantly between the full trial population and the genetic substudy. Whereas there were 14 (0.3%) reported deaths in the COLCORONA trial9, none of the participants to the genetic study died (P < 0.001), and there were 229 (5.1%) hospitalisation for COVID-19 in the trial compared to 58 (3.1%) in the genetic study (P < 0.001).

Table 1 Characteristics of the study population.
Full size table

GWAS for time to remission of COVID-19 symptoms

There were 1723 participants included in the GWAS for time to remission of COVID-19 symptoms. Of those, 1254 (72.8%) reported end of symptoms during the study period and 469 were censored at the date of their last telephone visit. We found two candidate genomic regions significantly associated with time to symptom remission located on chromosomes 5 and 9 (Fig. 1). Sex-stratified GWAS did not find additional associations. We found one significant association signal in the GWAS for time to COVID-19 symptom remission using both study arms at the 5p13.3 locus with the leading variant rs1173773 (P = 4.94 × 10–8; MAF 0.34). When conditioning on rs1173773, no additional genetic variants remained significant at P < 5 × 10–8 in the region and rs1173773 had the highest probability of being causal according to statistical and functional prioritization. The minor allele (C) was associated with a higher remission rate (hazard ratio (HR) = 1.25, 95% CI 1.15–1.35) as compared to the T allele, irrespective of treatment arm (interaction P = 0.18) (Table 2). According to survival modeling, 44%, 54% and 59% of patients with the TT, CT and CC genotypes respectively were predicted to have had symptom remission by day 15 of the study (Fig. 2, Table 3). Using an accelerated failure time model, we estimated the acceleration factor for the effect of the genetic variant on treatment remission which provides a more intuitive measure of the genetic effect than the hazard ratio. We calculated that the C allele contributes to symptom remission with an acceleration factor of 0.84 (95% CI 0.80–0.89) compared to the T allele (P = 5.23 × 10–9). Compared to patients with the rs1173773-TT genotype, patients with the CC genotype had symptom remission with an acceleration factor of 0.73 (95% CI 0.64–0.82; P = 3.89 × 10–7), and the CT genotype had an acceleration factor of 0.82 (95% CI 0.75–0.89; P = 6.30 × 10–6). The rs1173773 variant is located in intron 3 of the most abundant transcript of the natriuretic peptide receptor 3 gene (NPR3). We found no evidence supporting a regulatory role for this variant.

Figure 1
figure 1

Manhattan plots for the GWAS of time to remission of COVID-19 symptoms. (a) Using a Cox proportional hazards regression with 1723 subjects from the colchicine and placebo arms of the COLCORONA study, controlling for study arm, sex, age, and 10 principal components for genetic ancestry, with 6,392,715 genetic variants of minor allele frequency ≥ 5%. (b) Using a Cox proportional hazards regression with 851 subjects from the placebo arm of the COLCORONA study, controlling for sex, age, and 10 principal components for genetic ancestry, with 6,390,776 genetic variants of minor allele frequency ≥ 5%.

Full size image
Table 2 Genetic association results of the leading genetic variants identified in the GWAS for time to COVID-19 symptom remission.
Full size table
Figure 2
figure 2

Cumulative incidence curves for COVID-19 symptom remission.

Full size image
Table 3 Estimated cumulative incidence of COVID-19 symptom remission events by genotype groups.
Full size table

Results in placebo-treated participants

The most significant association in the GWAS for time to remission of COVID-19 symptoms in the 851 participants from the placebo arm was at 9q33.1 with the intergenic variant rs62575331 (P = 2.95 × 10–8, MAF = 0.12). When conditioning on rs62575331, no additional genetic variants remained significant at P < 5 × 10–8 in the region and rs62575331 had the highest probability of being causal according to statistical and functional prioritization. The minor allele (G) at the rs62575331 variant leading the 9q33.1 signal was associated with symptom remission (HR = 1.57, 95% CI 1.34–1.84). The association of the variant with symptom remission in the colchicine arm was not significant (P = 0.15) and the interaction between the genetic variant and colchicine treatment was significant (P = 1.19 × 10–5), highlighting a possible interaction with colchicine on symptoms duration that may be due to a shared biological pathway between the genetic variant and colchicine (Table 2). We calculated that 46%, 62% and 64% of those in the placebo-group with the CC, CG and GG genotypes respectively had symptom remission by day 15 of the study (Fig. 2, Table 3). In accelerated failure time modeling, compared to the C allele, the G allele contributed to symptom remission with an acceleration factor of 0.71 (95% CI 0.64–0.80; P = 4.83 × 10–9). Patients with the rs62575331-GG genotype had symptom remission with an acceleration factor of 0.52 (95% CI 0.36–0.75; P = 5.24 × 10–4), and patients with the CG genotype an acceleration factor of 0.71 (95% CI 0.61–0.82; P = 2.09 × 10–6) as compared to the CC genotype. The second-best GWAS signal in the GWAS with the placebo group did not reach statistical significance but was at locus 5p13.3, which was concordant to the region found by GWAS with participants from both study arms. The rs62575331 variant was not associated with the risk of hospitalisation for COVID-19 in the placebo arm (P = 0.87) or in the colchicine arm (P = 0.61) and had no interaction effect with colchicine on hospitalisation risk.

Sensitivity analyses

Sensitivity analyses were conducted with the two genetic variants identified by GWAS by adjusting for possible factors of COVID-19 disease severity, including BMI, obesity, diabetes, hypertension, smoking status, history of respiratory disease, and coronary artery disease. In each of the tested model, the addition of the factors did not impact the association of the genetic variant with time to symptom remission (all P values < 5 × 10–8). Our functional analysis of the GWAS results has highlighted that the two genetic variants were previously reported to be associated with height. In our study population, height was not a significant predictor of time to symptom remission (P = 0.67), the two genetic variants were not associated with height (P > 0.05), and the association of the two genetic variants with time to symptom remission was robust to further adjustment for height (P < 5 × 10–8).

GWAS for hospitalisation due to COVID-19

We conducted a GWAS for hospitalisation due to COVID-19 which occurred in 58 (3.1%) of the 1855 subjects included in the analysis. The study had limited power (Supplementary Methods) and none of the tested genetic variants passed the GWAS significance threshold (P < 5 × 10–8) (Supplementary Fig. 2).

Discussion

We report the results of a genome-wide association study of time to remission of COVID-19 symptoms in an outpatient population recently diagnosed with COVID-19 and who presented with at least one risk factor for COVID-19 complications. We found two genomic regions associated with symptom remission located at 5p13.3 and 9q33.1. The 5p13.3 locus spans the NPR3 gene, encoding a receptor for the binding of the natriuretic peptides which is involved in the clearance of natriuretic peptides, diuresis, blood pressure, and cardiometabolic diseases27. NPR3 is involved in the extended renin-angiotensin system which has been proposed as a possible mechanisms involved in the development of lung injury in COVID-1928. The rs1173773 C allele associated with remission of COVID-19 symptoms at the chromosome 5 locus has previously been associated with greater standing height29, and variants in NPR3 have also been associated with forced expiratory volume in 1 s (FEV1)29, and systolic blood pressure30.

The rs62575331-C allele associated with symptom remission at the leading variant at 9q33.1 was also previously associated with greater standing height29, similarly to variants in the pappalysin 1 gene (PAPPA) (P < 1 × 10–45)23,29,31 which is the nearest coding gene by transcriptional start site. The rs62575331 variant is an eQTL of a neighbouring lncRNA increasing the risk of venous abnormalities21,32. The 9q33.1 locus variant had a genetic effect on symptom remission which was limited to participants in the placebo arm, and it is possible that the pharmacological effects of colchicine share a common biological pathway with the effects of the rs62575331 variant. In the COLCORONA trial, colchicine was not found to modulate duration of COVID-19 symptoms. It is possible that the genetic variant at 9q33.1 may act in a biological pathway that is shared with that of colchicine, which acts on multiple inflammatory pathways including the NLRP3 inflammasome33,34. A recent report by Stella and colleagues35 on familial Mediterranean Fever (FMF) in the context of COVID-19 argues that there could be a finely regulated competition between the NLRP3 inflammasome and pyrin that is necessary to maintain protective inflammation levels. FMF is an inherited monogenic autoinflammatory disorder caused by excess activity of the pyrin protein and which is typically treated with colchicine to prevent fever attacks.

The genomic regions identified in this study are novel and do not overlap with previous genetic reports on host genetics of COVID-19 disease. So far, efforts to study the contribution of host genetics to COVID-19 have largely focused on the identification of genetic variants for risk of severe COVID-19 outcomes such as death or hospitalisation for COVID-19 and for genetic variants of risk of infection with the SARS-CoV-2 virus3,4,5,6. To the best of our knowledge, ours is the first to study the genetics of COVID-19 symptom remission in an outpatient population. The findings have yet to be replicated using an independent study sample. Concordant with our observations, a survey of the Center for Disease Control has reported that up to 35% of outpatient individuals still had symptoms 14 to 21 days after a COVID-19 diagnosis1, which included young adults without underlying chronic medical conditions. Symptom duration is an important predictor of disease severity which may involve immune responses, as supported by the higher T cell responses in mildly symptomatic individuals with a SARS-CoV2 infection as compared to asymptomatic individuals36. Some patients are experiencing prolonged multiorgan symptoms and complications beyond the initial period of acute infection and illness. A survey led by the UK Government's Office for National Statistics reported that one in five people who tested positive for COVID-19 had symptoms that lasted for 5 weeks or longer, and one in ten had symptoms that lasted for 12 weeks or longer2. There could be some shared pathophysiological mechanisms between time to COVID-19 symptom remission and post-acute COVID-19 syndrome37, and it may be relevant to assess the effect of the identified genetic variants in studies evaluating post-acute COVID-19 syndrome outcomes. Considering the important economic and health burden of both acute and post-acute COVID-19, genetics may provide insight toward the development of preventative interventions to lessen this burden.

Limitations

Our study had some limitations. Notably, we had limited power to study the more severe outcomes of death or hospitalisation for COVID-19 due in part to a healthy volunteer bias. The invitation and consent to the genetic substudy occurred after the randomisation visit, and very ill patients may have been less likely to agree to the additional interview or may have already died or been hospitalised by the time of recontact. We made all efforts to rapidly reach all participants and persisted recontact attempts even after the end of the 30-day follow-up period on treatment. Healthy volunteer bias is frequent in optional pharmacogenomic studies of clinical trials. To diminish this bias, a simultaneous consent process with the main study and the collection of genetic material at the randomisation visit are recommended. In practice, however, the additional consent for the genetic study adds time to the recruitment process and can be a deterrent to overall participation, particularly for a disease with acute onset such as COVID-19. Information on symptom types was not collected in this study. Because analyses were conducted with individuals of European genetic ancestry, validation of the genetic associations in other populations will be necessary. The inclusion and exclusion criteria of the COLCORONA trial may contribute to limiting the representativeness of the population. Importantly, the results have not yet been replicated in an independent population sample, and we cannot exclude the possibility that the results may be chance findings. As such, the results are not ready to guide clinical decisions.

Conclusion

This is the first study to report a GWAS for time to remission of symptoms in non-hospitalized patients with COVID-19 with at least one risk factor for a severe form of the disease. We found two genomic regions associated with symptom remission located at 5p13.3 and 9q33.1 in individuals diagnosed with COVID-19. The findings are novel and will need to be replicated in an independent sample.