Coronavirus infection, or COVID-19, caused by the SARS-CoV-2 virus, has resulted in one of the largest pandemics in human history (1). Epidemiological data on this disease emphasizes the remarkable heterogeneity in the course of the disease, ranging from completely asymptomatic cases to ICU hospitalizations with ventilator support and even fatal outcomes (2–5).

COVID-19 is associated with several comorbidities and non-genetic risk factors that increase the likelihood of developing the disease and experiencing severe progression, including respiratory failure. These risk factors include older age, male gender, and certain medical conditions such as cardiovascular disease, diabetes, obesity, chronic respiratory and kidney diseases, immunodeficiency, and neurological disorders (6–10).

The observed variability in susceptibility to and course of coronavirus infection suggests that both nongenetic risk factors and genetic variants may contribute to the clinical heterogeneity of COVID-19 in the population. Therefore, the COVID-19 Host Genetic Initiative (HGI) conducted a large study on a multiethnic sample of more than 49,500 COVID-19 patients from 46 studies in 19 countries (11). The study identified 13 genome-wide significant loci that are associated with SARS-CoV-2 infection or severe manifestations of COVID-19. Along with other genome-wide association studies (GWAS), several key genetic variants associated with severe COVID-19 outcomes were identified, which included those related to the function of the SARS-CoV-2 receptor (ACE2, ABO, TMPRSS2, and SLC6A20) (12–15) and immune response to the virus (HLA-region) (16, 17).

As epidemiological data and associated genetic variants for COVID-19 disease continue to accumulate, the combined assessment of disease severity for clinical implications remains an area of ongoing research. It is hypothesized that a combination of different gene variants may determine the severity of COVID-19 disease course. In other words, the complex polymorphic pathogenesis of coronavirus infection suggests a polygenic architecture.

Horowitz et al. (13) used the data from the COVID-19 HGI to establish a polygenic risk score (PRS) model. They demonstrated that individuals in the top 10% of the COVID-19 PRS among Europeans (n = 44.958) had a 1.38-fold increased risk of hospitalization and 1.58-fold increased risk of severe disease. Similarly, Farooqi et al. (18) reported a 1.57-fold higher risk of severe COVID-19 in high-risk patients compared to low-risk individuals. Another study by Crossfield et al. (19) involving 9,560 UK Biobank participants showed an adjusted odds ratio (OR) of 1.32, [95% confidence interval (CI): 1.11–1.58] for the highest PRS quintile compared with the lowest one. Nostaeva et al. (20) applied the same HGI statistics to a small cohort of Russian patients (1,085 participants, 347 individuals with severe COVID-19, and 738 with moderate or without disease) low-pass whole genome sequencing (LP-WGS). They found that more than one million genetic variants can be included in calculating polygenic scores to stratify patients by the risk of severe COVID-19. Individuals in the top 10% of the PRS distribution had over a two-fold increased risk of severe COVID-19 (odds ratio, OR: 2.2; 95% CI: 1.3–3.3, P-value = 0.0001).

Since coronavirus infection can develop rapidly in just over a week, identifying individuals at risk of developing severe COVID-19 through genetic variants may help identify targeting agents for investigating appropriate therapeutic interventions. Despite large-scale vaccination programs, optimal treatment selection remains a topical challenge. Therefore, the generation of a PRS model and patient profiling based on the risk of severe COVID-19 course can serve as a valuable tool for the healthcare industry.

In this research, we investigate the genetic factors that contribute to the severity of COVID-19 within the Eastern European population, which is underrepresented in many studies. By utilizing GWAS, we aim to identify genetic variants associated with COVID-19 and explore how these variants differ from those found in other populations. Additionally, we develop a PRS model to identify individuals at a high risk of experiencing severe COVID-19 outcomes.

To account for population stratification, we constructed a PCA plot with two principal components (Supplementary Figure S1) and performed clustering to detect minor populations and outliers. The largest cluster (number 1) was selected for further analysis and all other samples were considered outliers.

We performed a GWAS analysis on 7,124 individuals (with 53% being female individuals) from the selected cluster. Among them, 6,400 individuals were controls (i.e., they were aged over 40 years and had recovered from COVID-19 without experiencing a severe outcome), and the remaining 724 patients were cases (i.e., they experienced a severe course of COVID-19). The mean age was 64.77 years for cases and 49.53 years for controls. The characteristics of the final cohort are summarized in Table 1.

The GWAS analysis revealed several genome-significant loci (P-value < 5 × 10⁻⁸), (Table 2). Figure 1 shows the Manhattan plot with the GWAS results.

Particularly, the most genome-wide significant loci were identified in the leucine zipper transcription factor like 1 (LZTFL1) gene: rs35624553 (OR = 1.58, 95% CI: 1.34–1.85, P-value = 3.26 × 10⁻⁸) and rs10490770 (OR = 1.57, 95% CI: 0.082–1.34, P-value = 3.51 × 10⁻⁸). These findings were previously mentioned in the summary statistics of the COVID-19 HGI relevant to our phenotype. Three intergenic variants on chromosome 3 (rs17713054, rs13078854, and rs71325088) had the same association with severe forms of COVID-19 as observed with rs10490770.

Furthermore, we compared the effects of SNPs identified by the COVID-19 HGI (11), which are associated with critical illness (similar to our phenotype) and the current GWAS results. We could not compare the effect of rs912805253 as it is absent in hg19. The remaining SNPs exhibited poor approximation with HGI results; however, the intergenic variants highlighted in the HGI summary statistics (rs1819040 and rs74956615) demonstrated the effects that were comparable to our results (rs1819040: OR_hgi = 0.906, OR_genotek = 0.920, rs74956615: OR_hgi = 1.434, OR_genotek = 1.461). Figure 2 illustrates these results.

The heritability h2 coefficient was found to be 0.05 ± 0.0432.

Furthermore, we generated the PRS model to identify the group of individuals with a high risk of severe COVID-19 (Figure 3). In this group, the risk of experiencing a severe course of the disease is approximately twice that of the median risk [with OR = 2.25 (1.7, 2.95), P-value = 3.1 × 10⁻⁹]. The area under the curve (AUC) for the developed PRS is equal to 0.6 (0.58–0.62). The PRS was constructed using the grid model in LDPred2 and included 955,503 SNPs.

In this study, we identified several genome-wide significant loci associated with the severe forms of COVID-19. The loci include LZTFL1 and PRKCH-AS1/PRKCH genes, intergenic variants on chromosome 3, and a novel loci with intergenic variants on chromosome 6 with genome-wide significance (P-value = 5.92 × 10⁻⁸ to 6.15 × 10⁻⁸). In addition, several genetic variants in the LZTFL1 gene showed genome-wide significance.

The intron variant rs1989566, located within the PRKCH-AS1/PRKCH gene, was newly associated with the susceptibility and severity of COVID-19. The PKC family serves as a mediator for diverse signaling pathways and governs numerous crucial cellular functions, including proliferation, differentiation, and apoptosis. One of the PKC family members, protein kinase C eta protein (PKCη), which is encoded by the PRKCH gene, is a serine-threonine kinase. It is predominantly expressed in vascular endothelial cells and plays a role in the progression and exacerbation of atherosclerosis and subsequently stroke, as indicated by previous studies (29, 30).

Furthermore, the identified variant may potentially impact the expression of the PRKCH-AS1/PRKCH gene and subsequently disrupt endothelial function, which can cause an imbalance in hemostasis favoring a procoagulant state. This is characterized by impaired vasodilator release, increased release of vasoconstrictors, heightened microvasculature spastic reactions, increased leukocyte migration across the endothelium, and the initiation of localized inflammation. Prolonged exposure to factors that induce endothelial dysfunction can contribute to a pro-inflammatory and prothrombotic phenotype in endothelial cells (31). This exposure also results in a depletion of the pool of progenitor endothelial cells (32), ultimately limiting the capacity for restoring their normal phenotype and function.

In addition to regulating cell proliferation, differentiation, and cell death, PKCη is also expressed in the lung tissue, immune system, and proliferation pathways, for instance, in activating nuclear factor κB (NF-κB) signaling, which leads to anti-cancer drug resistance (33–35).

We confirmed the effect of rs10490770 on severe course of COVID-19 (OR = 1.57, P-value = 3.51 × 10⁻⁸) in the LZTFL1 gene. Previous studies (11, 36) have shown that rs10490770 increases the risks of all-cause mortality (HR = 1.4), severe respiratory failure (OR, 2.1), venous thromboembolism (OR = 1.7), and hepatic injury (OR = 1.5). For elderly people, it significantly increases the risk of mortality or severe respiratory failure by more than a double (OR, 2.7).

The LZTFL1 gene, expressed in the normal lung epithelium, is involved in protein transport to the cilia of the ciliated epithelium respiratory cells (37). The transcriptome analysis of lung biopsies from patients with COVID-19 showed the presence of signals associated with epithelial-mesenchymal transition of lung cells (EMT) or pulmonary fibrosis, which is regulated by LZTFL1, suggesting that this locus may serve as a potential therapeutic target (38). Other studies have also linked the LZTFL1 gene at the 3p21.31 locus to COVID-19 infection (39, 40).

Several limitations could affect our GWAS results, including relatively small number of controls, low trait heritability (h2 coefficient), and high polygenicity. Previous research has established that age is a significant risk factor for severe COVID-19, with older individuals being at a higher risk. In our study, the patient cases were sourced from hospitals with an average age of 64.77 years, while controls were selected from the database of a genetic testing company, Genotek, with an average age of ~35 years. We did not use a population control, as is commonly done in many studies. Instead, we defined controls as individuals aged 40 years or above who reported having COVID-19 but did not experience a severe course of the disease. The age threshold was set at 40 years to balance the age distribution and the sizes of our case and control cohorts.

Our primary aim was to evaluate the predictive ability of the PRS of COVID-19 severity in the Eastern European cohort. Therefore, we performed the GWAS analysis on a cohort of 7,124 individuals and constructed a polygenic risk model. Stratification of individuals based on PRS quantiles revealed that the high-risk category (top 10%) had twice the risk of severe course of COVID-19 compared to the median risk group. These results support findings from previous European cohort studies demonstrating similar associations of the PRS with severe COVID-19 (13, 18–20). Before implementation, the PRS must be validated in independent cohorts. Additional prospective studies can be beneficial to assess the clinical utility of the PRS in a practical setting. The potential application of the PRS could involve stratifying all tested individuals into high- and low-risk groups. In our study, we demonstrated the predictive power of the developed PRS. The selection of a threshold for the PRS that delineates the high-risk group should be determined based on factors such as the anticipated increase in hospitalization risk, mortality rates, and economic considerations.