Germline variants in UNC13D and AP3B1 are enriched in COVID-19 patients experiencing severe cytokine storms

Critically ill coronavirus disease 2019 (COVID-19) is characterized by severe cytokine storms, a hyperinflammatory condition intimately related to the development of fatal outcomes. Why some individuals seem particularly vulnerable to severe cytokine storms is still unknown. Primary immunodeficiency (PID)-related genes are inherited factors that dysregulate host inflammatory responses to infection, especially hemophagocytic lymphohistiocytosis (HLH)-related genes, established as contributors to the development of excessive cytokine storms. We analyzed the association between PID gene variants with severe cytokine storms in COVID-19. We conducted whole-exome sequencing in 233 hospitalized COVID-19 patients and identified four PID gene (UNC13D, AP3B1, RNF168, DHX58) variants were significantly enriched in COVID-19 patients experiencing severe cytokine storms. The total percentage of COVID-19 patients with variants in UNC13D or AP3B1, two typical HLH genes, was dramatically higher in high-level cytokine group than in low-level group (33.3 vs. 5.7%, P < 0.001). Germline variants in UNC13D and AP3B1 were associated with the development of severe cytokine storms, fatal outcomes in COVID-19. These findings advance the understanding of individual susceptibility to severe cytokine storms and help optimize the current management of COVID-19.


Introduction
Coronavirus disease 2019 (COVID-19) is triggered by severe acute respiratory coronavirus 2 (SARS-CoV-2), taking on various degrees of severity, ranging from mild to critical illness [1]. Critically ill or intensive care unit (ICU)admitted COVID-19 is characterized by severe cytokine storms, a hyperinflammatory condition intimately related to the development of lethal acute respiratory distress syndrome, multiple organ failure, and even death [2,3]. One question about the cytokine storms still lingers: why some individuals seem particularly vulnerable? Few relevant research studies are now underway exploring the underlying genetic mechanisms.

Patients and sample collection
The study population consisted of 233 participants with newly diagnosed COVID-19 between January 13, 2020 and March 5, 2020 at Tongji Hospital, Huazhong University of Science & Technology, Wuhan, China. Informed consent was obtained from each patient or patients' relatives. The detailed information was shown in Supplementary materials.

High-and low-level cytokine groups
The levels of serum cytokines were determined by Bio-Plex Pro Human Cytokines 48-Plex Screening (Bio-Rad Life Sciences, Hercules, CA, USA) on a Luminex FlEXMAP 3D system (Luminex, Austin, TX, USA) by following the manufacturer's protocols. The relative levels (fold change relative to healthy donors) of IP-10, IL-Ra, or MCP-3 were used to perform receiver operating characteristic (ROC) curve analysis to discriminate between ICU and non-ICU patients. ICU admission, as a status variable, was used in ROC analysis since COVID-19 patients admitted to ICU tended to experience severe cytokine storms.
Exome sequencing and rare, damaging variants filtering All coding variants were generated from whole-exome sequencing (WES) or targeted exome sequencing of DNA samples. We constructed a knowledge base of 241 genes involved in host antiviral immunity (Table S1). The rare variants (minor allele frequency (MAF) ≤ 0.01 in overall populations from public database 1000 G and gnomAD) were identified. The damaging variants meant having a functional impact of the proteins and were predicted by in silico analysis based on PolyPhen-2, SIFT, and CADD. Variants were classified as damaging by at least two software predicted as damaging. The analysis strategies in detail were shown in Supplementary materials.

Results
We performed WES on 233 hospitalized COVID-19 patients (mild, n = 58; severe, n = 133; critical, n = 42) with an attempt to understand the correlation between PID genes and severe cytokine storms in COVID-19. The median age was 60 years (interquartile range: 46-68 years). More clinical data are detailed in Table S2 and analysis strategies presented in Fig. S1.
The severity of cytokine storms in COVID-19 was rated in terms of three cytokines (IP-10, IL-1Ra, and MCP-3) as previously reported [11]. Patients were included in high-level cytokine group (n = 15) when all three cytokines were significantly increased (above the cutoff values), and the others were assigned to low-level cytokine group (n = 174), see Fig. 1A and Supplementary materials. There were no significant differences in the term of age, sex, comorbidity, duration of hospitalization between the high-and low-level cytokine groups, suggesting similar clinical baseline characteristics between the two groups (Table S3). We focused on PID genes by constructing a knowledge base of 241 antiviral immunity-related genes, which were mainly reported in PID diseases (see Supplementary materials and Table S1) [4]. We identified a total of 534 rare (MAF ≤ 0.01), damaging variants in 145 PID genes from 233 COVID-19 patients (listed in Table S4). To find out genes most relevant to the severe cytokine storms in COVID-19 patients, the variant burden of each identified genes was compared between high-and lowlevel cytokine groups. Variants in four genes were found to be significantly enriched in high-level cytokine group. UNC13D (26.7%, odds ratio (OR) = 6.7, 95% confidence interval (CI) = 1.8-25.1) was the most frequently mutated gene, followed by AP3B1 (13.3%, OR = 13.2, 95% CI = 1.7-101.7), RNF168 (13.3%, OR = 26.6, 95% CI = 2.3-313.3), and DHX58 (20.0%, OR = 6.0, 95% CI = 1.4-26.0). The variant burdens of these four genes were also strikingly higher in highlevel cytokine group than in Chinese general population (n = 301) from the 1000 Genomes project (Fig. 1B and Table S5). UNC13D (encoding Munc13-4) and AP3B1 (encoding AP3), two typical HLH genes, were among the major players in the development of HLH due to defects in T/NK-cell cytotoxic granule release pathway. Three heterozygous missense variants (Gly863Asp) variant in UNC13D is predictably affected the MHD2 domain by the HOPE software. In high-level cytokine group, one patient had two variants in the UNC13D gene, and four patients had a monoallelic variant in the UNC13D or AP3B1 gene. The total percentage of COVID-19 patients with variants in UNC13D or AP3B1 was dramatically higher in high-level cytokine group than in low-level group (33.3 vs. 5.7%, P < 0.001), in critical group than in non-critical group (16.7 vs. 5.2%, P < 0.01), in ICU than in non-ICU admission group (16.1 vs. 5.9%, P < 0.05), and in death group than in survival group (21.7 vs. 5.7%, P < 0.01) (Fig. 1C). To further validate our findings, we examined the variants in UNC13D and AP3B1 genes by targeted sequencing in an additional outpatient group including 51 asymptomatic SARS-CoV-2 infection individuals. Only 2.0 % (1 out of 51) of the individuals had the UNC13D variant, 3.9 % (2 out of 51) had AP3B1 variants (Table S6). The total percentage of patients with UNC13D and AP3B1 variant in the outpatient group was much lower than in high-level cytokine group (5.9 vs. 33.3%, P < 0.05) (Fig. S3).

Discussion
Our study found that three types of PID gene variants were enriched in COVID-19 patients experiencing severe cytokine storms. UNC13D and AP3B1 variants might be associated with the development of severe cytokine storms, critical illness, and fatal outcomes in COVID-19 patients, while the functional relevance of two other type genes (antibody production gene RNF168 and viral sensor Fig. 1 Germline variants were significantly enriched in COVID-19 patients with severe cytokine storms. A The red, blue, and green lines represent the receiver operating characteristic (ROC) curves of IP-10, IL-Ra, or MCP-3 to discriminate COVID-19 patients with ICU admission or without ICU admission, respectively. The areas under the curve (AUC) and the optimal cutoff values set by Youden's index are shown in the graph. B The bar chart shows the mutation percentage of the four genes in high-and low-level cytokine groups and in Chinese population from 1000 G. The statistically significant mutation burden was analyzed by the chi-square test and the Fisher's test. *P < 0.05; **P < 0.01. C The stacked bar shows the total percentage of patients harboring UNC13D or AP3B1 variants. The patients were grouped according to the severity of cytokine storms, COVID-19 severity, ICU admission, and fatal outcomes. The statistical significance between the groups was assessed by the chi-square test. *P < 0.05; **P < 0.01; ***P < 0.001. regulatory gene DHX58) warrants further investigation. A growing number of studies reported that heterozygous variants in HLH-related genes might be associated with HLH development [14][15][16]. Our study suggests that carrying a putatively deleterious variant in the UNC13D or AP3B1 gene might increase the risk of disease severity, possibly via interaction with other genetic mutations and/or nongenetic factors (e.g., exposure to SARS-CoV-2). Regarding the limitation of this study, we cannot differentiate functionally deleterious or neutral variants based on computational analysis. Except for the splicing site variant, all the other variants, including the in-frame variants with high CADD scores might be neutral when experimentally tested [17,18] and the burden test should be further performed on experimentally proven deleterious variants. Our findings will be helpful in advancing the understanding of individual susceptibility to severe cytokine storms and optimizing the current management of COVID-19. For vulnerable populations genetically subject to the severe cytokine storms, more aggressive immunomodulation therapies should be considered [19].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.