The TNFRSF13C H159Y Variant Is Associated with Severe COVID-19: A Retrospective Study of 500 Patients from Southern Italy

To identify host genetic determinants involved in humoral immunity and associated with the risk of developing severe COVID-19, we analyzed 500 SARS-CoV-2 positive subjects from Southern Italy. We examined the coding sequences of 10 common variable immunodeficiency-associated genes obtained by the whole-exome sequencing of 121 hospitalized patients. These 10 genes showed significant enrichment in predicted pathogenic point mutations in severe patients compared with the non-severe ones. Moreover, in the TNFRSF13C gene, the minor allele of the p.His159Tyr variant, which is known to increase NF-kB activation and B-cell production, was significantly more frequent in the 38 severe cases compared to both the 83 non-severe patients and the 375 asymptomatic subjects further genotyped. This finding identified a potential genetic risk factor of severe COVID-19 that not only may serve to unravel the mechanisms underlying the disease severity but, also, may contribute to build the rationale for individualized management based on B-cell therapy.

Advanced age, male gender, and chronic comorbidities can worsen the clinical disease course. However, the host genetic background also plays a role in disease variability. From this perspective, the determinants of different susceptibilities to SARS-CoV-2 mostly involve genes related to the initial stages of infection, as demonstrated for genetic variants of the ACE2 and TMPRSS2 genes [3][4][5][6]. Furthermore, common variants at diverse loci (3p21.31, 9q34.2, 19p13.3, 12q24.13, and 21q22.1) and inactivating rare mutations in genes belonging to the type I interferon pathway have been found to associate with severe COVID-19 [7][8][9].
Conversely, the determinants of variable COVID-19 severity mainly include components of the host immune response to the infection. Approximately 40 genes have been associated with a susceptibility to SARS-CoV-2 infection and 50% of them with disease severity. They are involved in multiple immunological pathways including either the innate or the adaptive immune responses [2,10]. Regarding an innate immunity, several studies highlighted the importance of the intact type I interferon pathway, as well as of immune signaling pathways, including Toll-like receptors [10]. Instead, for adaptive immunity, the impact of primary immunodeficiencies associated with antibody deficiencies on a COVID-19 outcome is not yet clarified. A recent survey on primary antibody deficiencies described a severe clinical course in patients with common variable immunodeficiency (CVID) compared to patients with pure agammaglobulinemia. This observation suggests a contribution of immune dysregulation, including dysfunctional B cells, to the severe inflammatory phenotype in CVID [11].
We herein investigated the host genetic factors associated with COVID-19 severity by analyzing a case series of 500 patients from Campania, a region of Southern Italy, where the pandemic caused 325,178 infections and 6816 deaths by May 2021. A targeted analysis of whole-exome sequencing (WES) data highlighted a mutational enrichment of CVIDassociated genes in severe patients compared to non-severe ones. In particular, severely affected subjects showed a recurrent rare variant, p.His159Tyr (H159Y), in the TNFRSF13C gene, encoding the B cell-activating factor receptor (BAFFR) [12]. These findings contribute to unraveling the mechanisms underlying the disease severity and open new avenues for targeted therapeutic approaches.

Collection and Clinical Stratification of the Patients
From March to September 2020, 500 patients affected by COVID-19 were enrolled from various public hospitals located in Campania, a region of South Italy.
In Campania, at the beginning of March 2020, the most common Pango lineage was the B.1 (74%, sample count: 66); whereas, at the end of September 2020, B.1.177 was the most common strain (38%, sample count: 26) (https://covid19dashboard.regeneron.com, accessed on 1 June 2021). Of note, none of these two strains are classified as Variantof-Interest or Variant-of-Concern by the Centers for Disease Control and Prevention (https://www.cdc.gov/coronavirus/2019-ncov/variants/variant-info.html, accessed on 1 June 2021). None of the enrolled patients were vaccinated at the time of the study.
Relevant features of the patients were summarized in Table 1. Among them, 121 (24.2%) hospitalized patients were subdivided into two age-and gender-matched subgroups: (i) 39 were severe cases that required mechanical ventilation (MV), and (ii) 82 were non-severe hospitalized patients that did not need MV. No differences in coexisting diseases between the two subgroups were observed. As a control group, we selected 379 individuals over 45 years of age that were positive for anti-SARS-CoV-2-IgG antibodies but did not show any symptoms. The antibody positivity was assessed by the Anti-SARS-CoV-2 Elecsys E2G 300 assay (Roche Diagnostics, Monza, Italy). A group of 283 healthy individuals (male 89.5%; mean age 45 ± 6.6 years) was also included in the study.

DNA Extraction, Quantification, and Library Preparation for Sequencing
Genomic DNA of COVID-19 patients was extracted from peripheral blood using a Maxwell ® RSC Blood DNA Kit (Promega, Madison, WI, USA), and DNA concentration and purity were evaluated using a NanoDrop™ 8000 spectrophotometer. For library preparation, we measured the DNA concentration with a Qubit ® DNA Assay Kit with a Qubit ® 2.0 fluorometer (Life Technologies by Thermo Fisher Scientific, Waltham, MA, USA), and 1.0 µg of DNA was used for library preparation. Genomic DNA was sonicated to generate 180-250-bp fragments and then end-polished, A-tailed, and ligated with the full-length adapter with further PCR amplification. PCR products were purified using the AMPure XP system (Beckman Coulter, Brea, CA, USA) and quantified using the Agilent high-sensitivity DNA assay from the Agilent TapeStation 4200 system (Agilent Technologies, Santa Clara, CA, USA). The whole exome was captured with Agilent SureSelect Human All Exome V6 (Agilent Technologies, Santa Clara, CA, USA), and the sequencing was performed on an Illumina NovaSeq 6000 System (Illumina, San Diego, CA, USA).

Bioinformatic Analysis of Sequencing Data
The paired-end (2 × 150 bp) sequencing returned an average of 42 million raw reads per sample. After removing the sequencing artifacts, we retained, on average, 98.35% of the raw reads. The percentages of the bases with quality scores above 20 and above 30 (Q20 and Q30) were 97.7% and 93.8%, respectively. The cleaned reads were aligned versus the reference genome (GRCh37/hg19; RefSeq assembly accession: GCF_000001405.13), and mapping BAM files were obtained with BWA-mem (version 0.7.17) and SAMTools (version 1.8). On average, 99.89% of reads were mapped, and 20.47% of duplicate reads were removed with Picard (version 2.18.9). We covered 99.63% of the target regions, and the average sequencing depth of the target was 149.06x. Ninety-five point eight percent of the target regions were covered with at least 20 reads, which were sufficient for reliable variant calls. The SNVs and small insertions and deletions (INDELs) were detected with GATK HaplotypeCaller, and the functional annotation of the variants was performed with ANNOVAR. We obtained 189,625 raw SNVs and 32,089 raw INDELs per sample. We excluded off-target variants (e.g., intergenic, intronic, etc.) and single-nucleotide polymorphisms (SNPs) with allele frequencies greater than 1% in non-Finnish European populations of the 1000 Genomes Project, ExAC (v3), and GnomAD (v2.1.1) databases. To remove possible false positives, we also discarded variants falling in the genomic duplicated regions. The set of exonic variants was then filtered to remove synonymous SNVs.

SNP Genotyping
The asymptomatic patients and the healthy controls were genotyped by the TaqMan ® SNP Genotyping Assay (C__89561246_10) for the single-nucleotide variant (SNV) rs61756766 (Applied Biosystems by Thermo Fisher Scientific, Waltham, MA, USA).
Allele and genotype frequencies were compared using the chi-square, Fisher's exact, or Armitage tests. A two-sided p < 0.05 was considered statistically significant.

Mutational Enrichment Analysis
The overall features of the 500 COVID-19 patients enrolled in the study are summarized in Table 1.
To identify the SNVs significantly enriched in severe COVID-19 patients compared to non-severe ones, we performed a targeted analysis of the WES data using an in silico panel composed of 10 genes associated with CVID, as derived from the Human Phenotype Ontology database (Table 2).

Discussion
As for other viral infections, we can broadly classify the immune response to SARS-CoV-2 into two phases. At the early stage, a specific adaptive immune response is required to eradicate the virus and to prevent disease progression to severe stages. At the severe stage, the cytokine release syndrome and hyperactivation of the immune system are the main causes of life-threatening respiratory disorders leading to severe disease and, in some cases, to death [13].
Genetic differences are well-known to contribute to individual variations in the immune response to pathogens [2]. A few reports about the role of humoral immunity in COVID-19 have been published so far. A recent study in patients with primary antibody deficiencies and concurrent COVID-19 showed that patients with CVID (i.e., the presence of dysfunctional B lymphocytes) presented with a severe form of the disease requiring MV. Conversely, patients with agammaglobulinemia (i.e., the absence of B lymphocytes) presented mild symptoms. Thus, the different clinical course of COVID-19 in these two types of patients suggests a possible role of B lymphocytes in SARS-CoV-2-induced inflammation [11].
To understand the role played by the host genetic determinants involved in humoral immunity, we performed a WES analysis focused on CVID-associated genes by comparing SNV enrichment between severe and non-severe patients. Interestingly, we identified a rare variant in TNFRSF13C that was recurrent among severely affected patients. The TNFRSF13C gene encodes the BAFFR that is activated by BAFF, a TNF family member that supports the survival of B cells [12]. The c.475C > T transition results in the H159Y amino acid change within the highly conserved cytoplasmic domain of BAFFR [14]. This variant has already been associated with immune disorders, such as CVID, autoimmune lymphoproliferative syndrome-like disease, multiple sclerosis, Good's syndrome, and Sjogren's syndrome [15][16][17][18][19][20]. Additionally, it has been also associated with an increased risk of chronic lymphocytic leukemia and lymphoma [14,21].
Despite that the variant was not found to be associated with a change in BAFFR mRNA or protein expression [22], it has been proven to be a gain-of-function variant. Indeed, it increases ligand-independent BAFFR signaling, leading to enhanced NF-kB1 and NF-kB2 activation [14]. It is well-known in the role of NF-kB signaling in the induction of proinflammatory genes, including those encoding cytokines and chemokines [23]. Of note, it has been recently demonstrated that SARS-CoV-2 infection also induces an inflammatory phenotype in alveolar type 2 cells mediated by the hyperactivation of NF-kB signaling [24]. Thus, it is conceivable to hypothesize that the severe clinical course observed in COVID-19 patients carrying the rs61756766 AG genotype could be due to enhanced NF-kB signaling activation.
One limitation of this study might be the relatively low number of cases. However, it should be considered that our cohort of cases was derived from a careful selection to define homogeneous groups of patients. Furthermore, all the studied individuals were of Italian origin and residents in the same geographical area, which limited the problem relative to the population structure. Replication studies are needed to further validate our findings. Additionally, we herein focused on the H159Y variant, since it is the most frequently detected rare SNV in severe patients and since it happens to have strong biological support from the literature. Nevertheless, other non-tested variants could be associated with severe COVID-19. Thus, we encourage other research groups to carry out further genetic studies considering our preliminary results.
Our findings are relevant not only from the perspective of the pathogenesis of this condition but, also, for therapeutic implications. Indeed, the role of B cells in determining inflammatory disorders is also supported by the evidence of the beneficial effects of B-cell depletion in the management of these diseases [11,20]. Thus, the identification of the humoral immunity risk factors associated with a worse COVID-19 outcome could provide a rationale for individualized management based on B-cell therapy.  Informed Consent Statement: Informed consent was obtained from all the subjects involved in the study.

Data Availability Statement:
The data that supported the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.