This study included 241 hospitalized COVID-19 patients recruited from Huoshenshan hospital at Wuhan city, Hubei province, China between 11 January 2020 and 11 March 2020. COVID-19 was diagnosed based on chest computed tomography (CT) manifestations and/or reverse transcription-polymerase chain reaction (RT-PCR) following the criteria of the New Coronavirus Pneumonia Prevention and Control Program (5th edition). In this study, the mild COVID-19 patients were those with no obvious clinical symptoms or with fever, respiratory symptoms, and radiological evidence of pneumonia. The severe COVID-19 patients were those having at least one of the following conditions: respiratory distress, respiratory rate ≤30 beats/minute; mean oxygen saturation ≤93% in a resting state; arterial blood oxygen partial pressure/oxygen concentration ≤300 mm·Hg; respiratory failure and requiring mechanical ventilation; shock; and admission to intensive care unit (ICU) with other organ function failure. Classifications of COVID-19 severity were taken as the worst classification during the patient’s hospital stay.

The clinical characteristics of the patients were extracted from the electronic medical records. We collected three broad classes of characteristics: 1) demographic variables (age, sex, and ethnicity); 2) symptoms (fever and diarrhea); and 3) comorbid conditions (hypertension, diabetes, cardiac disease, chronic bronchitis, chronic liver disease, chronic obstructive pulmonary disease, cerebrovascular disease, and cancer).

The Ethics Committee of Huoshenshan Hospital approved the study (HSSLL036). Given the urgency of the COVID-19 pandemic, the need for informed consent was waived by the ethnics boards of the hospital.

We collected 30 unique variants associated with susceptibility or severity of SARS-CoV/SARS-CoV-2 infection from 21 papers (up to 29 November 2020), referred as “literature dataset”. In addition, we downloaded COVID19-hg GWAS meta-analysis results (release 4) produced by COVID-19 host genetics initiative (HGI) from https://www.covid19hg.org/results/r4/, and 1420 unique variants were selected that meet one of the following four requirements: 1) variants with p < 1E-5 in “A1_ALL” group, that is, phenotype of very severe respiratory confirmed covid vs. not hospitalized covid; 2) variants with p < 1E-5 in “B1_ALL” group, that is, phenotype of hospitalized covid vs. not hospitalized covid; 3) variants with p < 5E-7 in any one of group; 4) variants with p < 1E-5 in any three of groups except “A1_ALL” and “B1_ALL”. This was referred as “HGI dataset”. Two variants were overlapped in the two datasets, resulting in a total of 1448 unique variants collected.

IFNs play a central role in innate immunity against virus infection (Zhang et al., 2020a; Bastard et al., 2020) and cell receptors for virus are key determinant for viral entry (Hoffmann et al., 2020a; Shang et al., 2020). Therefore, We collected 19 IFN-pathway genes previously implicated in SARS-CoV/SARS-CoV-2 susceptibility and severity (Hamano et al., 2005; He et al., 2006; Ching et al., 2010; Zhang et al., 2020a; Zhang et al., 2020b) as well as 3 human cell receptors/co-receptors for SARS-CoV-2 (Hoffmann et al., 2020a; Hoffmann et al., 2020b; Zhou et al., 2020b; Shang et al., 2020). A total of 25 potential enhancers for these genes were obtained from GeneHancer database (GH score > 1, Gene Association > 100) (Fishilevich et al., 2017).

RNA probes were designed to cover these variants and gene enhancers. All the 25 enhancers and 1238 unique candidate variants were covered, including 29 variants in the literature dataset and 1211 variants in the HGI dataset (Supplementary Tables S1–S3).

Peripheral whole blood samples were collected from all participants. Genomic DNAs were extracted from 1 ml of peripheral whole blood, according to the manufacturer’s instructions (QIAamp DNA blood kits). The quality of the isolated genomic DNA was verified by the following two methods: 1) the DNA degradation and contamination were monitored in 1% agarose gels; and 2) the DNA concentration was measured using a Qubit DNA Assay Kit and a Qubit 3.0 Fluorometer (Life Technologies).

The targeted capture sequencing was conducted by iGeneTech Bioscience Corporation (Beijing, China). Briefly, Human genomic DNA was sheared to 150–200 bp by Bioruptor Pico (Diagenode). Then end repair, dA-tailing, and adapter ligation were performed. The ligation product was cleaned up and size-selected by using Beckman Ampure XP Beads (Beckman). The purified ligated product was amplified by using PCR. Then the library was in solution hybrid with biotinylated RNA probes, captured with Dynabeads MyOne Streptavidin T1 (Invitrogen), and amplified with PCR. The library was quantified by Qubit and fragment-size measured by Agilent 2100 Bioanalyzer system before high-throughput sequenced by NovaSeq.

Raw reads were firstly quality trimmed with Trimmomatic (Bolger et al., 2014). Clean reads were then aligned to the human reference genome (hg38) using BWA algorithm (Li and Durbin, 2009). PCR duplicates were removed using samtools (Li et al., 2009), and GATK software (McKenna et al., 2010) was used to call SNPs and indels. The detected variants were finally saved as VCF files. Data of autosomal biallelic variants for Han Chinese of 1000 Genomes Project (Auton et al., 2015) were downloaded from https://www.internationalgenome.org/. Genetic association analysis was conducted using PLINK 1.9 software (Chang et al., 2015).

p-values comparing demographics severe and mild disease groups were calculated by means of χ2 or Fisher exact test as appropriate, except for p-value for age which was calculated using student’s t test. Minor allele frequency (MAF) of variants was compared between case and control groups using Fisher exact test. Logistic regression analysis was also used to compute the contributing variants to severity with adjusting age, sex and comorbidities. Statistic power was calculated using G*Power 3 software (Faul et al., 2007). p < 0.05 was considered statistically significant.

A total of 241 COVID-19 patients were recruited, with 96 mild cases and 145 severe cases. Demographic and phenotypic data are shown in Table 1. Comparing mild patients with severe ones, the median age increased from 58 to 67 years (p < 0.0001). In addition, we observed a greater percentage of severe patients with comorbidities than mild patients (p = 0.03), specifically with diabetes (p = 0.04) and cerebrovascular disease (p = 0.03). This is in accordance with previous report that older COVID-19 patients and those with comorbidities were more likely to be severe in disease (Zhou et al., 2020a; Richardson et al., 2020). Previous report also indicated that severe patients had a greater percentage of male cases (Richardson et al., 2020). However, no obvious sex difference was found between mild and severe patients in our data (p = 0.17).

The experimental design was illustrated in Figure 1. To validate previously reported variants associated with COVID-19 susceptibility and severity, we collected 30 unique variants from 21 papers (referred as “literature dataset”) and selected 1420 unique variants from COVID-19 host genetics initiative (HGI) release 4 (Initiative, 2020) (referred as “HGI dataset”). Additionally, to identify regulatory variants associated with COVID-19 susceptibility and severity in the enhancers of previously reported associated genes, we collected 19 IFN-pathway genes previously implicated in coronavirus susceptibility and severity (Hamano et al., 2005; He et al., 2006; Ching et al., 2010; Zhang et al., 2020a; Zhang et al., 2020b) as well as 3 human cell receptors/co-receptors for SARS-CoV-2 (Hoffmann et al., 2020a; Hoffmann et al., 2020b; Zhou et al., 2020b; Shang et al., 2020) and obtained 25 potential enhancers for these genes from GeneHancer database (Fishilevich et al., 2017) (referred as “enhancer dataset”).

A panel of RNA probes was designed to capture these variant and enhancers, resulting in 1238 unique variants and all the 25 potential enhancers covered (Supplementary Tables S1–S3). The 1238 unique variants included 29 variants in the literature dataset and 1211 variants in the HGI dataset, with two variants overlapped between the two datasets. Variants that fail to design probes might be due to GC content, repetitive sequence, dimer or secondary structure of probes.

The probe panel was used to targeted capture sequence the 241 COVID-19 patients. A median of 8.5 million raw reads were obtained for each sample. After filtering, a median of 8.3 million reads was kept as clean reads, with a median of average insert size of 173 bp. The reads were mapped to hg38 genome. The median mapping rate was 99.56%, with a median duplication rate of 27.47% (Supplementary Figure S1). Specifically, the median number of target mapped reads was 4.5 million, with a median target coverage rate of 99.57% and median target mean depth of 1472× (Supplementary Figure S2).

Out of the 1238 candidate variants in the panel, a total of 1006 variants in the panel were identified in the 241 COVID-19 patients, including 26 variants in the literature dataset and 982 variants in the HGI dataset.

Comparing 96 mild and 145 severe COVID-19 patients, seven variants were confirmed to be significantly different in minor allele frequency (MAF) (p < 0.05, Supplementary Table S2), corresponding to four independent loci. The four lead variants were shown in Table 2. Specifically, rs1799964 in the promoter of inflammation-related gene TNF was found to be associated with COVID-19 severity. CC genotype of this variant has been reported to be associated with femoral head necrosis after SARS-CoV infection (Wang et al., 2008). In addition, rs2224986 and rs13062942 were still associated with COVID-19 severity after adjusting age, sex and comorbidities on regression analysis. Both variants had a significant difference between hospitalized and non-hospitalized patients in the HGI dataset.

When comparing all COVID-19 patients with ancestry-matched general population from 1000 Genomes Project, we identified 39 significant variants (p < 0.05, Supplementary Table S5), corresponding 10 independent loci. The independent lead variants were shown in Table 3. Specifically, multiple genetic variants in the ABO gene locus and 3p21.31 region have been validated. In addition, the missense variant rs429358 in exon of APOE gene, which has been reported to be associated with COVID-19-positive in the UK Biobank data (Kuo et al., 2020), was confirmed in Chinese population. As for the confirmed variant rs1886814 in the intron of FOXP4-AS1 (forkhead box P4 antisense RNA 1), recent trans-ethnic genome-wide association study of severe COVID-19 that incorporated Chinese population and HGI results also revealed a significant variant nearby, rs1853837 (Wu et al., 2021), which is LD with rs1886814 (1000 Genomes Project CHB, r ² = 0.68).

When comparing severe COVID-19 patients with ancestry-matched general population from 1000 Genomes Project, we identified 23 significant variants (p < 0.05, Supplementary Table S6), corresponding 7 independent loci. The independent lead variants were shown in Table 4. Again, multiple genetic variants in the ABO gene locus and 3p21.31 region, rs429358 in exon of APOE gene and rs1886814 in the intron of FOXP4-AS1 have been validated. What’s more, rs9975538 in the intron of gene IFNAR2 was also confirmed. IFNAR2 gene, along with IFNAR1 encodes type I interferon receptor, thus the rs9975538 might affect type I interferon pathway.

In total, 49 unique variants were confirmed to be associated with susceptibility or severity of COVID-19 (p < 0.05), corresponding 18 independent loci. It is worth noting that, for the confirmed variants, the phenotypes of the cases and controls were highly consistent between our study and previous reports, and the confirmed variants identified between mild and severe patients were quite different from those identified between patients and general population, suggesting the genetic basis of susceptibility and severity of SARS-CoV-2 infection might be quite different.

Previous candidate gene study of genetic association with COVID-19 mainly focused on the exon region of important genes (Zhang et al., 2020a; Novelli et al., 2020). However, regulatory regions, particularly enhancers, play a vital role in gene expression and may affect disease susceptibility and severity when misfunction (Claringbould and Zaugg, 2021). To investigate the role of enhancer variants of key genes in COVID-19, we also included 25 enhancers of 16 IFN-pathway genes and 3 SARS-CoV-2 receptor/co-receptor genes in the probe panel. These enhancers were predicted by GeneHancer database. Genetic association analysis was conducted between mild and severe COVID-19 patients, between all COVID-19 patients and ancestry-matched general population from 1000 Genomes Project, or between severe COVID-19 patients and ancestry-matched general population from 1000 Genomes Project.

In total, we identified 67 variants in the enhancer region of the panel that were associated with susceptibility or severity of COVID-19, relating to 12 enhancers of 11 genes (p < 0.05, Supplementary Tables S7–S9). Further annotation by RegulomeDB database filtered out five potential regulatory variants for broad-spectrum antiviral gene IFITM3 and one for SARS-CoV-2 co-receptor gene TMPRSS2 (probability score > 0.8, Table 5). Specifically, among the six variants, GTEx revealed two variants affecting the expression of IFITM3, that is, T allele of rs11246060 in the enhancer of IFITM3, which protected COVID-19 patients from severe outcomes when comparing mild and severe patients (p = 2.38E-2, OR = 0.39, 95% CI = 0.17–0.89), was associated with increased expression of IFITM3 in whole blood (p = 1.72E-10), while T allele of rs28655829 in the enhancer of IFITM3, which increased the risk of severity when comparing mild and severe patients (p = 2.38E-2, OR = 0.39, 95% CI = 0.17–0.89), was associated with decreased expression of IFITM3 in cultured fibroblasts (p = 7.55E-8). This indicated that these variants might confer genetic risk or protection by affecting the gene expression and highlighted the importance of IFITM3 gene in the defense of SARS-CoV-2 infection.