DOCK2 is involved in the host genetics and biology of severe COVID-19

Abstract

Identifying the host genetic factors underlying severe COVID-19 is an emerging challenge^1-5. Here we conducted a genome-wide association study (GWAS) involving 2,393 cases of COVID-19 in a cohort of Japanese individuals collected during the initial waves of the pandemic, with 3,289 unaffected controls. We identified a variant on chromosome 5 at 5q35 (rs60200309-A), close to the dedicator of cytokinesis 2 gene (DOCK2), which was associated with severe COVID-19 in patients less than 65 years of age. This risk allele was prevalent in East Asian individuals but rare in Europeans, highlighting the value of genome-wide association studies in non-European populations. RNA-sequencing analysis of 473 bulk peripheral blood samples identified decreased expression of DOCK2 associated with the risk allele in these younger patients. DOCK2 expression was suppressed in patients with severe cases of COVID-19. Single-cell RNA-sequencing analysis (n = 61 individuals) identified cell-type-specific downregulation of DOCK2 and a COVID-19-specific decreasing effect of the risk allele on DOCK2 expression in non-classical monocytes. Immunohistochemistry of lung specimens from patients with severe COVID-19 pneumonia showed suppressed DOCK2 expression. Moreover, inhibition of DOCK2 function with CPYPP increased the severity of pneumonia in a Syrian hamster model of SARS-CoV-2 infection, characterized by weight loss, lung oedema, enhanced viral loads, impaired macrophage recruitment and dysregulated type I interferon responses. We conclude that DOCK2 has an important role in the host immune response to SARS-CoV-2 infection and the development of severe COVID-19, and could be further explored as a potential biomarker and/or therapeutic target.

Fig. 1. GWAS in a Japanese population stratified by COVID-19 severity and age.

a, Forest plots of the risk of COVID-19-associated variants in a Japanese population. Error bars indicate the 95% confidence interval. b, Manhattan plot of the GWAS on severe COVID-19 in young patients (those less than 65 years of age) (440 cases and 2,377 controls). Uncorrected P values from the GWAS analysis are shown. The dotted line represents the genome-wide significance threshold of P < 5.0 × 10⁻⁸. Manhattan and quantile–quantile plots of all GWAS results are presented in Extended Data Fig. 2. MT, mitochondrial. c, Regional association plot at the DOCK2 locus. Dots represent SNPs coloured according to linkage disequilibrium (r²) with the lead SNP of rs60200309. FAM196B is also known as INSYN2B. d, Allele frequency spectra of the rs60200309-A allele in the 1000 Genomes Project Phase3v5 database.

Fig. 2. Cell-type- and tissue-specific expression of DOCK2 and its downregulation in severe COVID-19.

a, eQTL effect of the COVID-19 risk variant (rs60200309) on DOCK2 expression levels using bulk RNA-seq of peripheral blood. The risk allele (rs60200309-A) decreases DOCK2 levels in patients with COVID-19 aged below 65 years. TPM, transcripts per kilobase million. b,c, Differential expression analysis of DOCK2 with varying COVID-19 severity. DOCK2 expression levels were quantified by qPCR and normalized to GAPDH expression. b, Comparison between severe and non-severe COVID-19 cases. c, Comparison between most severe, severe, mild and asymptomatic cases of COVID-19. d–k, scRNA-seq in PBMCs from individuals with severe COVID-19 (n = 30) and healthy controls (n = 31). d, Uniform manifold approximation and projection (UMAP) visualization of all 394,526 cells. e, Projection of DOCK2 gene expression. Innate immune cell clusters are outlined with a red dashed line. f, Percentage of DOCK2-expressing cells and DOCK2 expression levels. g, Expression change with severe COVID-19 in six major cell types. h, Visualization and annotation of the innate immune cell clusters. i–k, DOCK2 expression and expression changes with severe COVID-19 in the innate immune cell clusters. i, Projection of DOCK2 gene expression. j, Percentage of DOCK2-expressing cells and DOCK2 expression levels. k, Expression change with severe COVID-19 in five cell types. l, COVID-19 context-specific decreasing eQTL effect of the DOCK2 risk variant in non-classical monocytes. m,n, Immunohistochemical analysis of DOCK2. Lung and hilar lymph nodes were obtained from patients with COVID-19 pneumonia (m) or controls without COVID-19 or pneumonia (n), and stained with anti-DOCK2 polyclonal antibody. Results for all samples are shown in Extended Data Fig. 9. Scale bars, 0.25 mm. In a–c,l, boxes denote the interquartile range (IQR) and the median is shown as horizontal bars; whiskers extend to 1.5 times the IQR; outliers are shown as individual points in a–c and all samples are shown as individual points in l. Uncorrected P values are shown in (a–c,g,k,l). cDC, conventional dendritic cells; cMono, classical monocytes; intMono, intermediate monocytes; Mono, monocytes; ncMono, non-classical monocytes; NK, natural killer cells; Pro T, proliferative T cells; T_reg, T regulatory cells.

Fig. 3. In vivo suppression of DOCK2 in a Syrian hamster model of SARS-CoV-2 infection.

a, Changes in body weight of hamsters infected with SARS-CoV-2. b, Representative images of lungs collected after euthanizing the hamsters at 11 dpi. c, Lung weight changes after infection. The number of samples (n) is indicated. d, Representative lung histopathology and immunohistochemistry of the infected hamsters at 6 dpi. Outlined areas are expanded to the right of each image. Right, lung tissue was stained with the anti-CD68 mouse monoclonal antibody to highlight alveolar macrophages. e, SARS-CoV-2 viral loads in the organs of the infected hamsters at 3 and 6 dpi. f, Lung cytokine expression assays of the infected animals. Ip-10 is also known as CXCL10. In (a) and (c), the error bars represent standard error of the mean, and P values were determined with two-sided Welch’s t-test; *P < 0.05; **P < 0.01; ***P < 0.001. In (e) and (f), boxes denote the IQR, and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR, and all animals are shown as individual points. P values were determined with two-sided Wilcoxon rank sum test.

Extended Data Fig. 1. Japan COVID-19 Task Force.

Japan COVID-19 Task Force is a nation-wide consortium to overcome COVID-19 pandemic in Japan, which was established in early 2020. Japan COVID-19 Task Force consists of > 100 hospitals (red dots) led by core academic institutes (blue labels), and collected DNA, RNA, and plasma from the COVID-19 cases along with detailed clinical information. The figure was originally created using sf and ggplot2 R packages based on Global Map Japan version 2.1 Vector data (Geospatial Information Authority of Japan).

Extended Data Fig. 2. A principal component analysis plot of the GWAS participants and Manhattan and quantile-quantile plots of the GWAS.

(a, b) A principal component analysis (PCA) plot of the GWAS participants (COVID-19 cases and controls) along with and without International HapMap populations (a and b, respectively). (c) Manhattan plots and quantile-quantile plots of the Japanese GWAS of COVID-19. Uncorrected P values from GWAS analysis are shown. Dotted lines represent the genome-wide significance threshold of P < 5.0 × 10⁻⁸.

Extended Data Fig. 3. Regional association plots of the HLA imputation analysis.

Regional association plots of the HLA imputation analysis results. Dots represent SNPs and HLA variants with colors according to the legend. Uncorrected P values from HLA imputation analysis are shown. Dotted lines represent the genome-wide significance threshold of P < 5.0 × 10⁻⁸. HLA genes with the most significant associations in each of the case-control phenotypes are indicated.

Extended Data Fig. 4. ABO blood type associations with COVID-19 in Japanese and cross-population Mendelian randomization analysis of the COVID-19 GWAS.

(a) Odds ratios of the ABO blood types in the Japanese population are indicated. Dots represent the odds ratios and bars represent the 95 % confidence intervals. P values are uncorrected. Detailed results are presented in Supplementary Table 5. (b) Forest plots of the Mendelian randomization (MR) analysis results of causal inference on the COVID-19 GWAS in Japanese (left panel) and Europeans (right panel). Since effect sizes (= beta) of MR are not scalable among phenotypes and populations, normalized beta is indicated. For each phenotype and population, the standard error for the COVID-19 GWAS with the largest sample size (i.e., “COVID-19 vs control” for Japanese and “Self-reported COVID-19 vs control (C2)” for Europeans) was set to be 0.1. Dots represent the effect size normalized beta estimates and bars represent the 95 % confidence intervals. P values are uncorrected. The abbreviations of the exposure phenotypes and the detailed MR results are given in Supplementary Table 6 and Supplementary Table 7. BMI; body mass index, T2D; type 2 diabetes, CPD; cigarettes per day, CAD; cardiovascular disease, SBP; systolic blood pressure, DBP; diastolic blood pressure, eGFR; estimated glomerular filtration rate, UA; serum uric acids, RA; rheumatoid arthritis, SLE; systemic lupus erythematosus.

Extended Data Fig. 5. Effect size comparisons of the COVID-19 risk loci between the discovery GWAS and the replication study.

Co-plots of the odds ratios and 95% confidence intervals between the discovery GWAS cohort and replication cohort. To focus on the differences in the cases collected in different pandemic waves (initial waves for GWAS and latter waves for the replication), same controls as GWAS were currently used for the cases in the replication. A regression coefficient was estimated based on logarithm of odds ratios. Dots represent the odds ratios and bars represent the 95 % confidence intervals.

Extended Data Fig. 6. Colocalization analysis of the GWAS and eQTL signals at the DOCK2 locus.

Regional colocalization plots of the GWAS signals (severe and younger COVID-19 cases vs controls) and the eQTL signals on DOCK2 expression in the COVID-19 patients at the DOCK2 locus. CLPP; colocalization posterior probability. The eQTL effects of the variants around DOCK2 region are given in Supplementary Table 10.

Extended Data Fig. 7. Cell type definition and gene ontology enrichment analysis of DOCK2 co-expression gene module in the PBMC single cell analysis.

(a) Violin plots showing the expression distribution of selected canonical cell markers in the 12 clusters of PBMC. The rows represent selected marker genes and the columns represent clusters with the same color as in Fig. 2d. (b) Violin plots showing the expression distribution of selected canonical cell markers in the 5 clusters of innate immune cell clusters, shown in the same color as in Fig. 2h. (c) Tile plot showing percentage concordance between the manually annotated 12 clusters and Azimuth annotation. (d) The top 25 enriched biological processes by gene ontology (GO) analysis of DOCK2 co-expression gene module identified by weighted gene co-expression network analysis (WGCNA) in the non-classical monocytes of COVID-19 patients, where DOCK2 showed the highest cell type-specific expression profile. The color of the dots represents the adjusted P values.

Extended Data Fig. 8. Biological impacts of DOCK2 downregulation in primary cells and DOCK2 knockdown and Interferon-α production assay in THP-1 Blue ISG cells.

(a) The impact of DOCK2 downregulation on interferon-α (IFN-α) production ability in pDC. Sorted pDC were stimulated with CpG and/or CPYPP. Data shows means ± s.e.m. (n = 3 per group). Differences of IFN-α production ability between the groups were evaluated using two-sided paired t-test. (b) The impact of DOCK2 downregulation on chemotaxis in CD3⁺ T cells. CD3⁺ T cells were stimulated with CXCL12 or CXCL12 + CPYPP (n = 19 per group). Differences of chemotaxis between the groups were evaluated using two-sided paired t-test. (c, d) Knockdown of DOCK2 by CRISPR system was confirmed by western blotting (c) and qRT-PCR. (d) Semi-quantitative staining density measure was determined using ImageJ (NIH). Data shows means ± s.e.m. (n = 3 per group). Data are compared to control group. P values were determined with One-way ANOVA followed by Dunnett’s post hoc test. (e, f) Activity ratio of SEAP reporter to no treatment group. Reporter was activated by 50 ng/ml LPS (e) or 50 μg/ml polyIC (f). Data shows means ± s.e.m. (n = 3 per group). Data are compared to control group. P values were determined with One-way ANOVA followed by Dunnett’s post hoc test.

Extended Data Fig. 9. Immunohistochemical analysis for DOCK2.

Lung and hilar lymph nodes were obtained from autopsied cadaver (Sample 1–3, 6, 7) or surgical specimen (Sample 4, 5), and stained by anti-DOCK2 polyclonal antibody. Sample 1–3; COVID-19 pneumonia. Sample 4-5; control. Sample 6; non-COVID-19 severe pneumonia. Sample 7; non-COVID-19 mild pneumonia.

Extended Data Fig. 10. In vivo suppression of DOCK2 in a Syrian hamster model with SARS-CoV-2 infection.

(a) Schematic timeline of the experimental procedure. (b) Changes in weight of uninfected animals. The error bars represent standard error of the mean. (c) Changes in weight of each of the infected animals, corresponding to Fig. 3a. Three CPYPP-administrated animals reaching humane endpoint were euthanized at dpi 7 and 9, lowering survival rate to 77% (=10/13), while survival of vehicle-administrated animals was 100% (=12/12). The animals were administered with CPYPP (red), or vehicle (blue). (d) Histopathological examination of the lungs of infected hamsters. Syrian hamsters were inoculated with SARS-CoV-2 with CPYPP or Vehicle. Syrian hamsters infected with CPYPP or Vehicle were euthanized on dpi 3, 6, and 11 for pathological examinations (n = 3). Shown are pathological findings in the lungs of hamsters infected with the virus on dpi 3, 6, and 11 (hematoxylin and eosin staining). Middle and Right show enlarged views of the area circled in black in Left. (Scale bars, 2.5 mm [Left], 1.0 mm [Middle], and 0.25 mm [Right].) (e) Immunohistochemistry for alveolar macrophages. Shown are immunohistochemical findings in the lungs of hamsters infected with the virus on dpi 6 (n = 3 per group). Lung tissue was stained with the anti-CD68 mouse monoclonal antibody. (Scale bars, 0.25 mm.) (f) Pathological severity scores in infected hamsters. To evaluate comprehensive histological changes, lung tissue sections were scored based on (d) pathological changes. Scores were determined based on the percentage of inflammation area of the maximum cut surface collected from each animal in each group by using the following scoring system: 0, no pathological change; 1, affected area (≤10%); 2, affected area (<50%, > 10%); 3, affected area (<90%, ≥50%); 4, (≥90%) an additional point was added when pulmonary edema and/or alveolar hemorrhage was observed. The total score is shown for individual animals. Blue dot and red dot indicate +Vehicle and +CPYPP, respectively.