The study included 199 hospitalized COVID-19 patients from Central Russia. The Ethical Review Committee of Kursk State Medical University approved the study protocol (protocol No. 1 from January 11, 2022), and all participants provided written informed consent. The inclusion criteria for the study required participants to have self-declared Russian ancestry and to have been born in Central Russia. Table 1 provides the baseline and clinical characteristics of the study cohort.

The patients were enrolled in the study during the COVID-19 pandemic from 2020 to 2022 at the intensive care units (ICU) of Kursk Regional Hospital №6 and Kursk Regional Tuberculosis Dispensary. All patients had a PCR-confirmed COVID-19 diagnosis. All patients underwent a complete clinical examination. Qualified doctors established the diagnosis of concomitant diseases based on a full range of clinical and laboratory-instrumental research methods. In our case–control study, all patients with severe COVID-19 were divided into subgroups depending on the presence or absence of T2DM or obesity.

In accordance with the WHO guidelines, low fruit and vegetable consumption was defined as consuming less than 400 g per day. Adequate consumption of fresh vegetables and fruits was defined as consuming 400 g or more, equivalent to 3–4 servings per day, excluding starchy tubers, such as potatoes. Insufficient physical activity was characterized by engaging in less than 180 minutes per week of moderate to vigorous physical activities. This encompassed various forms of exercise, including leisure activities such as walking and running, and fitness club exercises, such as treadmill running, aerobics, or resistance training.

When selecting genes and polymorphisms, we used data from the largest GWAS meta-analysis studying severe COVID-19, published in 2023 [7]. Our study selected loci with a p-level significance of $\le$1 $\times$ 10${}^{-19}$ established in this meta-analysis. Then, SNPs with a minor allele frequency 0.05 were excluded from the analysis, as well as loci for which our laboratory was unable to develop probes for genotyping using the TaqMan-based-PCR method for methodological reasons (low CG composition, presence of GC clamps, runs of identical nucleotides). In total, 10 SNPs were included in the genotyping: rs143334143 CCHCR1, rs111837807 CCHCR1, rs17078346 SLC6A20–LZTFL1, rs17713054 SLC6A20–LZTFL1, rs7949972 ELF5, rs61882275 ELF5, rs12585036 ATP11A, rs67579710 THBS3, THBS3-AS1, rs12610495 DPP9, rs9636867 IFNAR2.

Genotyping of the SNPs was performed using an allele-specific probe-based polymerase chain reaction (PCR) according to the protocols designed in the Laboratory of Genomic Research at the Research Institute for Genetic and Molecular Epidemiology of Kursk State Medical University. The Primer3 software (https://primer3.ut.ee/) was used for primer design. A real-time PCR procedure was performed in a 25 µL reaction solution containing 1.5 units of Hot Start Taq DNA polymerase (Biolabmix, Novosibirsk, Russia), approximately 10 ng of DNA, and the following concentrations of reagents: 0.25 µM of each primer; 0.1 µM of each probe; 250 µM of each dNTP; 3 mM MgCl${}_2$ for rs7949972, 3.5 mM MgCl${}_2$ for rs61882275, 2 mM MgCl${}_2$ for rs12610495, and 2.5 mM MgCl${}_2$ for the remaining SNPs; 1 × PCR buffer (67 mM Tris–HCl, pH 8.8, 16.6 mM (NH${}_4$)${}_2$SO${}_4$, 0.01% Tween-20). The PCR procedure comprised an initial denaturation for 10 minutes at 95 °C, followed by 39 cycles of 92 °C for 30 s and 57 °C, 59 °C, 60 °C, 61 °C, 62 °C, 63 °C, 65 °C, 66 °C for 1 min (for rs12610495, rs17078346, rs17713054, rs111837807, rs9636867, rs143334143 and rs7949972, rs12585036 and rs61882275, rs67579710, respectively). To assure quality control, 10% of the DNA samples were genotyped twice, blinded to the case–control status. Over 99% of the data were concordant.

STATISTICA software (v13.3, TIBCO Software Inc., Palo Alto, CA, USA) was utilized for statistical processing. The normality of the distribution for quantitative data was assessed using the Shapiro–Wilk test. Most quantitative parameters exhibited deviations from the normal distribution, so they were presented as the median (Me) along with the first and third quartiles (Q1 and Q3). The Kruskal–Wallis test was used to compare quantitative variables among three independent groups. Following that, groups were contrasted pairwise using the Mann–Whitney test, which was also performed to compare quantitative variables among two independent groups. For categorical variables, differences in statistical significance were evaluated using Pearson’s chi-squared test with Yates’s correction for continuity.

The compliance of genotype distributions with Hardy–Weinberg equilibrium was evaluated using Fisher’s exact test. The study groups’ genotype frequencies and associations with disease risk were analyzed using the SNPStats software (https://www.snpstats.net/start.htm). The log-additive model was considered for the genotype association analysis. Associations within the patients with/without obesity or T2DM were adjusted for covariates. Variables that revealed differences in the general clinical characteristics of the study groups were used as covariates.

The model-based multifactor dimensionality reduction (MB-MDR) analysis tested two-, three-, and four-level genotype combinations (G–G) and genotype–environment combinations (G–E). The empirical p-value (p${}_{\mathrm{\text{perm}}}$) was estimated using a permutation test for each model. Models with p${}_{\mathrm{\text{perm}}}$ 0.05 were considered statistically significant. All calculations were adjusted for covariates. Statistical analysis was carried out using the R software environment. Models (on average, 3–4 models for each level) with the highest Wald statistics and the lowest p-level significance were included in the study. Additionally, using the MB-MDR method, individual combinations of genotypes associated with the studied phenotypes were established (p 0.05). Calculations were performed in the MB-MDR program for the R software environment (Version 3.6.3) (https://www.r-project.org/) [24].

The MDR method also analyzed the most significant G–G and G–E models. The analysis was implemented in the MDR program (v.3.0.2) (http://sourceforge.net/projects/mdr). The MDR method was used to assess the mechanisms of interactions (synergy, antagonism, additive interactions (independent effects)) and the strength of interactions (the contribution of individual genes/environmental factors to the entropy of a trait and the contribution of interactions). The results of the MDR analysis were visualized as a graph.

The functional effects of SNPs were examined using bioinformatics resources, the methodologies and functionalities of which were extensively detailed in our previous studies [25, 26, 27, 28, 29].

• The bioinformatic tool GTExportal (http://www.gtexportal.org/) was used to analyze the expression levels of the studied genes in whole blood, blood vessels, adipose tissue, and pancreas, as well as the expression quantitative trait loci (eQTLs) [30].

• To examine SNPs binding to quantitative expression trait loci (eQTL) in peripheral blood, the eQTLGen resource available at https://www.eqtlgen.org/ was employed [31].

• HaploReg (v4.2), a bioinformatics tool available at https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php, was utilized to assess the associations between SNPs and specific histone modifications indicative of promoters and enhancers. These modifications included acetylation of lysine residues at positions 27 and 9 of the histone H3 protein, as well as mono-methylation at position 4 (H3K4me1) and tri-methylation at position 4 (H3K4me3) of the histone H3 protein. Additionally, the tool was applied to investigate the positioning of SNPs in DNase hypersensitive regions, regulatory motif sites, and locations binding to regulatory proteins [32].

• The atSNP Function Prediction online tool (http://atsnp.biostat.wisc.edu/search) was used to evaluate the impact of SNPs on the gene affinity to transcription factors (TFs) depending on the carriage of the reference/alternative alleles [33]. TFs were included based on the degree of influence of SNPs on the interaction of TFs with DNA and calculated based on a positional weight matrix.

• Using the Gene Ontology online tool (http://geneontology.org/), it was feasible to analyze the joint involvement of TFs linked to the reference/SNP alleles in overrepresented biological processes directly related to the pathogenesis of obesity and T2DM [34]. Biological functions controlled by transcription factors associated with SNPs were used as functional groups.

• Bioinformatic analyses were conducted utilizing the resources provided by the Common Metabolic Diseases Knowledge Portal (https://hugeamp.org/) and the Type 2 Diabetes Knowledge Portal (https://t2d.hugeamp.org/). These portals integrate and analyze data from the largest consortiums dedicated to studying metabolic diseases, facilitating comprehensive investigations into genetic associations, intermediate phenotypes, and risk factors related to obesity and T2DM.

The analysis of associations of genotypes in the entire group (Table 2) revealed a link between rs17713054 SLC6A20–LZTFL1 and rs7949972 ELF5 and the increased risk of obesity. Specifically, the rs17713054 SLC6A20–LZTFL1 demonstrated an odds ratio (OR) = 2.34, 95% confidence interval (CI) = 1.24–4.4, p = 0.007 (risk allele A), while the rs7949972 ELF5 exhibited an OR = 1.79, 95% CI = 1.11–2.91, p = 0.015 (risk allele T).

Subsequently, in the entire group of patients (Table 3), rs9636867 IFNAR2 was associated with a higher risk of T2DM (risk allele G, OR = 8.28, 95% CI = 1.69–40.64, p = 0.027).

Using the MB-MDR approach, the six most significant gene–gene interaction patterns associated with obesity in patients with severe COVID-19 were identified: Three two-locus, two three-locus, and one four-locus (p${}_{\mathrm{\text{perm}}}$ $\le$ 0.05) (Table 4). In total, the best models of G–G interactions included five polymorphic loci: rs7949972, rs17713054, rs61882275, rs12585036, and rs143334143; all of them participated in the formation of two or more of the most significant G–G interactions. In the next step, we analyzed the mechanisms of interactions between these genetic variants using the multifactor dimensionality reduction (MDR) method (Fig. 2).

First, MDR showed that the genetic variants included in the best G–G models are characterized predominantly by moderate and pronounced synergism, except for SNP rs17713054, which showed additive (independent) effects (in interaction with rs61882275 and rs143334143); in addition, additive effects were observed when analyzing the interactions of rs61882275 and rs7949972. Secondly, the effects of intergenic interactions (trait entropy indicators: 0.03%–4.60%) were comparable to the mono-effects of SNPs (entropy indicators: 0.21%–3.45%). Thirdly, the following combinations of genotypes of polymorphic gene variants have the strongest correlations with obesity in patients with severe COVID-19: rs7949972 ELF5 C/C$\times$rs17713054 SLC6A20–LZTFL1 G/G (Beta = –0.25103; p = 0.002), rs618822758 ELF5 G/G $\times$ rs17713054 SLC6A20–LZTFL1 G/G (Beta = –0.24582; p = 0.002), rs12585036 ATP11A C/T $\times$ rs17713054 SLC6A20–LZTFL1 G/G (Beta = –0.22951; p = 0.006), rs7949972 ELF5 C/C $\times$ rs12585036 ATP11A C/T $\times$ rs17713054 SLC6A20–LZTFL1 G/G (Beta = –0.31162; p = 0.012), rs618822758 ELF5 A/G $\times$ rs12585036 ATP11A C/T $\times$ rs143334143 CCHCR1 G/G (Beta = –0.38917; p = 0.005), rs618822758 ELF5 A/G $\times$ rs12585036 ATP11A C/T $\times$ rs17713054 SLC6A20–LZTFL1 G/G $\times$ rs143334143 CCHCR1 G/G (Beta = –0.53689; p = 0.001) (Supplementary Table 1).

Using the MB-MDR approach, the six most significant gene–gene interaction patterns associated with type 2 diabetes mellitus in patients with severe COVID-19 were identified: Three three-locus and three four-locus models (p${}_{\mathrm{\text{perm}}}$ $\le$ 0.05). None of the two-locus models showed a p-level significance at p $\mathrm{\le }$ 0.05 after performing the permutation test (Table 5). In total, the best models of G–G interactions included eight polymorphic loci, six of which, rs7949972, rs61882275, rs12585036, rs143334143, rs67579710, rs12610495, were involved in the formation of two or more of the most significant G–G interactions. In the next step, we analyzed the mechanisms of interactions between these genetic variants using the multifactor dimensionality reduction (MDR) method (Fig. 3).

First, MDR showed that genetic variants included in two or more best G–G models are characterized mainly by moderate and pronounced synergism, except for SNP rs143334143, which showed additive (independent) effects (in interaction with rs67579710 and rs12610495), as well as rs67579710 (in interaction with rs143334143 and rs61882275). Secondly, the effects of intergenic interactions (maximum entropy indicators of a trait: 4.15%) exceeded the mono-effects of SNPs (maximum entropy indicators: 1.56%). Thirdly, the following combinations of genotypes of polymorphic gene variants have the strongest correlations with T2DM in patients with severe COVID-19: rs9636867 IFNAR2 G/G $\times$ rs67579710 THBS3, THBS3-AS1 G/G $\times$ rs12585036 ATP11A C/T (Beta = –0.25861; p = 0.003), rs7949972 ELF5 T/T $\times$ rs67579710 THBS3, THBS3-AS1 G/G $\times$ rs12585036 ATP11A C/C (Beta = 0.31938; p = 0.011), rs7949972 ELF5 C/C $\times$ rs12610495 DPP9 G/G $\times$ rs143334143 CCHCR1 G/G (Beta = 0.400926; p = 0.004), rs67579710 THBS3, THBS3-AS1 G/G $\times$ rs618822758 ELF5 G/G $\times$ rs12610495 DPP9 G/G $\times$ rs11183780 CCHCR1 T/T (Beta = 0.49877; p = 0.001), rs7949972 ELF5 C/C $\times$ rs67579710 THBS3, THBS3-AS1 G/G $\times$ rs12610495 DPP9 G/G $\times$ rs143334143 CCHCR1 G/G (Beta = 0.498775; p = 0.001), rs67579710 THBS3, THBS3-AS1 G/G $\times$ rs618822758 ELF5 G/G $\times$ rs12610495 DPP9 G/G $\times$ rs143334143 CCHCR1 G/G (Beta = 0.50013; p = 0.001) (Supplementary Table 2).

Using the MB-MDR approach, the nine most significant models of gene–environment interactions associated with obesity in patients with severe COVID-19 were identified: Two two-level models, four three-order models, and three four-level models (p${}_{\mathrm{\text{perm}}}$ $\le$ 0.01) (Table 6). In total, the best G–E interaction models included physical inactivity interacting with SNPs. It is noteworthy that five genetic variants were found to be involved in two or more of the best G–E models: rs17713054, rs61882275, rs12585036, rs143334143, and rs7949972. In the next step, we analyzed the interactions between these genetic variants and physical activity level, fruit and vegetable intake, and smoking using the multifactor dimensionality reduction (MDR) method (Fig. 4).

First, MDR found that the environmental risk factor level of physical activity exhibited moderate antagonism in interaction with smoking and fresh fruit and vegetable consumption. Secondly, environmental risk factors are characterized by both moderate/pronounced synergism and additive (independent) effects in interaction with genetic variants included in two or more of the most significant models of gene–environment interactions. Thirdly, the mono-effects of environmental factors (contribution to the entropy of the trait is 1.62%–2.46%) are comparable to the mono-effects of SNPs (contribution to the entropy of the trait is 0.21%–3.45%). Fourth, the effects of gene–environment interactions (contribution to entropy: 0.09%–3.42%) are comparable to the mono-effects of environmental factors and SNPs included in the most significant G–E interactions (contribution to trait entropy: 0.21%–3.45%). Fifthly, the following gene–environmental interactions have the strongest correlations with obesity in patients with severe COVID-19: PHYS 1 $\times$ rs7949972 ELF5 T/C (Beta = 0.1712; p = 0.04), VEGET 1 $\times$ rs17713054 SLC6A20–LZTFL1 A/G (Beta = 0.22898; p = 0.01), SMOKE 0 $\times$ rs618822758 ELF5 A/A $\times$ rs12610495 DPP9 G/A (Beta = 0.44601; p = 0.008), SMOKE 0 $\times$ rs7949972 ELF5 C/C $\times$ rs17713054 SLC6A20–LZTFL1 G/G (Beta = –0.21285; p = 0.021), SMOKE 0 $\times$ rs618822758 ELF5 A/A $\times$ rs143334143 CCHCR1 G/G (Beta = 0.2946; p = 0.011), SMOKE 0 $\times$ rs618822758 ELF5 G/G $\times$ rs17713054 SLC6A20–LZTFL1 G/G (Beta = –0.20831; p = 0.024), VEGET 1 $\times$ rs618822758 ELF5 A/G $\times$ rs12585036 ATP11A C/T $\times$ rs143334143 CCHCR1 A/G (Beta = 0.48492; p = 0.017), SMOKE 0 $\times$ rs9636867 IFNAR2 A/G $\times$ rs17078346 SLC6A20–LZTFL1 A/A $\times$ rs12585036 ATP11A C/T (Beta = –0.36649; p = 0.003), VEGET 1 $\times$ rs618822758 ELF5 A/G $\times$ rs12585036 ATP11A C/C $\times$ rs11183780 CCHCR1 T/T (Beta = 0.26523; p = 0.02) (Supplementary Table 3).

Using the MB-MDR approach, the three most significant models of gene–environment interactions associated with the development of T2DM in patients with severe COVID-19 were identified, and after a permutation test, a statistically significant level of p $\mathrm{\le }$ 0.05 was established only for the three three-level models (Table 7). In total, the best G–E interaction models included physical inactivity interacting with five SNPs. In the next step, we analyzed the interactions between these genetic variants and physical activity levels using the multifactor dimensionality reduction (MDR) method (Fig. 5).

Firstly, MDR revealed that the mono-effect of physical inactivity (0.93% of trait entropy) is comparable to the mono-effects of SNPs of genes included in the most significant gene–environment interactions (0.41%–1.57% of entropy). Secondly, the effects of gene–environment interactions (maximum entropy values: 6.78%) exceed the mono-effects of physical inactivity (0.93%) and the mono-effects of SNPs (0.41%–1.57% entropy). Thirdly, physical inactivity is characterized by moderate and pronounced synergism in interaction with all genetic variants included in the best G–E models, except for interactions with rs67579710 and rs61882275, in which physical inactivity was characterized by independent (additive) effects. Fourthly, the effects of interactions of genetic variants included in the most significant G–E models (maximum entropy indicators: 4.31%) are comparable to the effects of genotype–environment interactions (maximum entropy indicators: 6.78%). Fifthly, the following gene–environment interactions have the strongest correlations with T2DM in patients with severe COVID-19: PHYS 1 $\times$ rs67579710 THBS3, THBS3-AS1 A/G $\times$ rs12585036 ATP11A C/T (Beta = 0.79917; p = 2.68 $\times$ 10${}^{-5}$) PHYS 1 $\times$ rs7949972 ELF5 T/T $\times$ rs12610495 DPP9 A/A (Beta = 0.509078; p = 2.44 $\times$ 10${}^{-2}$), PHYS 1 $\times$ rs618822758 ELF5 A/G $\times$ rs12610495 DPP9 A/A (Beta = 0.43616; p = 0.000482) (Supplementary Table 4).