Introduction
Metabolic diseases, such as type 2 diabetes mellitus (T2DM) and obesity, are important comorbidities in COVID-19 and are associated with higher risks of hospitalization, critical care admission, and death.1,2 A common feature of severe COVID-19, diabetes mellitus, and obesity is dysregulation of the immune system and a persistent inflammatory response.3 In COVID-19, immune dysfunction can manifest as acute and excessive inflammation, known as cytokine storm.4 This response is driven in part by monocyte hyperactivation, which can lead to excessive release of inflammatory cytokines.5 Recent studies confirm that elevated glucose levels can promote macrophage and monocyte activation through hypoxia-inducible factor-1α (HIF-1α), which upregulates the transcription of proinflammatory cytokines.6 Moreover, hyperglycemia increases angiotensin-converting enzyme 2 (ACE2) expression, which may increase monocyte susceptibility to viral entry and support SARS-CoV-2 replication.6,7 High glucose levels also impair several components of the innate immune system, such as natural killer cells’ activity, neutrophil migration, and phagocytosis,8,9 potentially weakening antimicrobial host defense.
Along with T2DM, obesity is characterized by a chronic, low-grade inflammatory state.10 Hypertrophic adipocytes recruit monocytes and promote their polarization toward an inflammatory phenotype, thereby promoting increased production of proinflammatory cytokines.3 Furthermore, adipokines greatly contribute to the hyperinflammatory state in obese individuals.11 For instance, leptin triggers proliferation of monocytes and macrophages, thereby increasing levels of inflammatory mediators.12,13 Because leptin concentrations correlate with adipose tissue volume,14 hypertrophy and proliferation of adipocytes may lead to excessive inflammation. Moreover, adipocytes express ACE2; therefore, adipose tissue may contribute to susceptibility to SARS-CoV-2 entry.15 Notably, adipocytes can also serve as reservoirs for several viral agents, such as adenovirus, influenza A virus, and cytomegalovirus, raising the possibility that SARS-CoV-2 could persist in adipose tissue and subsequently spread to other susceptible tissues.1
Prolonged inflammation, oxidative stress, hyperglycemia, and hyperlipidemia are significant stress factors for the cell. Heat shock proteins (HSPs) are a pivotal class of chaperones that protect cells against various exogenous and endogenous stressors.16-18 HSPs maintain proteostasis under stress conditions and regulate inflammatory pathways and may thereby limit cytokine-induced apoptosis.19 Furthermore, HSPs perform immunomodulatory functions. They are involved in antigen presentation, innate and adaptive immune responses,20 and cytokine production.21,22 HSPs also influence inflammatory pathways via ubiquitination, promoting or preventing proteolytic degradation of their members.23
The objective of this study was to analyze associations between polymorphic variants in genes encoding HSPs and T2DM or obesity among patients with severe COVID-19.
Materials and methods
This study was conducted as a case-control analysis within a cohort of patients with severe or critical COVID-19 and is reported in accordance with the STROBE guidelines. It was a case-control genetic association study within this consecutively enrolled COVID-19 cohort.
The 195 patients with severe or critical COVID-19 included in the present study were consecutively enrolled in 2022, after ethical approval, from the intensive care units (ICUs) of Kursk Regional Hospital No. 6 and Kursk Regional Tuberculosis Dispensary. The study population was derived from the same severe/critical COVID-19 cohort described in our previous publications.20,24,25 Four participants from the previously reported 199-patient cohort were excluded because their T2DM status was unavailable. COVID-19 diagnosis was confirmed by real-time polymerase chain reaction (PCR) for SARS-CoV-2 upon admission. Severe and critical COVID-19 were classified according to the Russian Ministry of Health guidelines in effect at the time of enrollment. Severe COVID-19 was characterized by high-grade fever, marked dyspnea (respiratory rate > 30 breaths per minute), SpO2 < 93%, a PaO₂/FiO₂ ratio < 300 mmHg, arterial lactate > 2 mmol/L, and radiographic findings typical of viral pneumonia. A critical course was diagnosed in the presence of COVID-19-specific complications: acute respiratory distress syndrome (ARDS), acute respiratory failure requiring invasive mechanical ventilation, multiple organ failure, and septic shock. Baseline and clinical characteristics of the study population are listed in Supplementary Table. 1 and 2.
No formal pre-study power calculation was performed because of limited regional case volume and restricted biobank sample availability. Participants were consecutively recruited patients with severe or critical COVID-19 from a two-center ICU cohort. Eligibility with respect to ancestry was based on self-reported Russian descent and birthplace; genetic ancestry or population structure was not formally evaluated using genetic principal components.
Among the study participants, 39 patients had T2DM, 101 had obesity.24,25 T2DM and obesity were analyzed as internal case-control outcomes. For the T2DM analysis, subjects with T2DM were considered cases and those without T2DM were considered controls. For the obesity analysis, subjects with obesity were considered cases and those without obesity were considered controls. The present study evaluates a different panel of HSP-related single nucleotide polymorphisms (SNPs), and the HSP SNP–T2DM/obesity association analyses reported here have not been published previously.
The study design is summarized in Figure 1.
The study was conducted in accordance with the Declaration of Helsinki, and received ethical approval from the Ethical Review Committee of Kursk State Medical University (Protocol No. 1, dated January 11, 2022). All participants provided written informed consent prior to enrollment.
T2DM diagnostic criteria included glycated hemoglobin A1c of 6.5% or higher, fasting plasma glucose of 7.0 mmol/L or higher, 2-h plasma glucose during a 75-g oral glucose tolerance test of 11.1 mmol/L or higher, or, in patients with classic symptoms of hyperglycemia or hyperglycemic crisis, random plasma glucose of 11.1 mmol/L or higher.26 It is important to note that T2DM status was determined based on a pre-existing diagnosis documented in the medical record and was not assigned during the acute COVID-19 admission, which allowed established T2DM to be distinguished from transient stress-related or glucocorticoid-related hyperglycemia. For the purposes of this study, subjects with a body mass index (BMI) of 30 kg/m² or higher were classified as having obesity.27 A low intake of fruits and vegetables was defined as daily consumption of fewer than 400 g.28 Insufficient physical activity was defined as less than 150 min of physical activity per week.29
Genotyping was performed in the Laboratory of Genomic Research at the Research Institute for Genetic and Molecular Epidemiology of Kursk State Medical University. Genomic DNA was extracted from venous blood samples using phenol/chloroform extraction and ethanol precipitation. The blood samples were collected into EDTA-coated tubes and stored at −20°C until processing. The quality and concentration of extracted DNA were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).
SNP genotyping was performed using TaqMan-based quantitative PCR according to protocols established at the Research Institute for Genetic and Molecular Epidemiology. Details of primer design, the reaction mixture, and PCR conditions were described previously.20 To ensure quality control, 10% of DNA samples were genotyped in duplicate, blinded to case-control status.
The selection of SNPs was based on their biological relevance, predicted functional potential, and technical feasibility. The analysis focused on genes encoding proteins in the HSP40, HSP70, HSP90, HSF, and BAG families. For these groups, we applied the following criteria: a minor allele frequency (MAF) of 0.05 or higher in Europeans, strong predicted functional impact according to bioinformatic analyses, and practical methodological constraints. SNPs were excluded based on the following criteria: substantial difficulties for primer and probe design (for instance, extreme GC content), close proximity to other variants or location within repetitive genomic regions, deviation from Hardy-Weinberg equilibrium (HWE) (p-value HWE <0.05), MAF below 0.05 in the study sample, low duplicate concordance rate. For HWE testing Pearson’s chi-squared criteria were applied. HWE analysis was a predefined quality control procedure conducted on the entire cohort (195 people) before the T2DM and obesity association analyses. Supplementary Table. 3 presents quality-control metrics and SNP-specific regression sample sizes for the genotyped HSP SNPs, including MAF, HWE results, duplicate concordance rates, and the number of participants included in each regression analysis. SNP rs6457452 HSPA1B deviated from Hardy-Weinberg equilibrium (p-value HWE <0.05) and was therefore excluded from further analysis. The remaining SNPs met the above criteria and were included in the further analysis.
Statistical evaluations, including quality control, genotype frequency assessments, and association analyses with T2DM and obesity among patients with severe COVID-19, were performed using PLINK version 1.9.0-b.7.7. All studied SNPs were tested for each outcome; T2DM and obesity were analyzed as separate outcomes. We applied an adaptive permutation test (Pperm), which reports pointwise empirical P values (EMP1). Statistical significance was defined at a threshold of Pperm < 0.05. This approach allows nonparametric evaluation of significance, thereby reducing the risk of inflated type I error. No family-wise or other across-SNP multiple-testing correction was applied. In addition, a log-additive model was utilized to estimate genotype associations. Because environmental and endogenous factors may modify genetic associations, we conducted adjusted analyses. For the T2DM analysis, models were adjusted for sex, smoking status, BMI, and fresh fruit/vegetable intake; for the obesity analysis, models were adjusted for sex, smoking status, T2DM status, and fresh fruit/vegetable intake. For functional analyses, we used a wide range of bioinformatic resources. To interpret functional effects, we included values with adjusted false discovery rate (FDR) < 0.05.
Statistical analyses were conducted using STATISTICA software (version 13.3). Data normality was evaluated using the Shapiro-Wilk test. Because most quantitative variables deviated from a normal distribution, results are reported as medians and quartiles (Me [Q1; Q3]). Group comparisons involving three independent samples were performed using the Kruskal-Wallis test, followed by pairwise Mann-Whitney U tests. Additionally, the Mann-Whitney test was applied to compare quantitative data between two independent groups.
To explore possible molecular mechanisms through which the studied SNPs may be related to severe COVID-19, we used a wide range of open bioinformatics resources. Functional annotation was performed only for polymorphisms that showed statistically significant associations with the outcomes. Only statistically significant annotations meeting the resources’ internal adjustment criteria were retained; no additional adjustment was applied to these external annotation results. The bioinformatic analysis was conducted entirely in silico and should not be interpreted as direct evidence of gene expression or protein-level change.
Bioinformatic resources included:
The GTExPortal (http://www.gtexportal.org/ ) and eQTLGen (https://www.eqtlgen.org/ ) were used to evaluate gene-expression patterns in tissues relevant to COVID-19, T2DM, and obesity, including whole blood, adipose tissue, pancreas, liver, skeletal muscle, blood vessels, and lungs and to analyze expression quantitative trait loci (eQTLs).30-32
The atSNP affinity assay (http://atsnp.biostat.wisc.edu/search) was utilized to evaluate how SNPs influence the binding affinity of DNA sequences to transcription factors (TFs). The atSNP performs affinity tests of motif matches with the SNP or the reference sequence and SNP-led changes in motif matches.32
The Gene Ontology web platform (http://geneontology.org/) offers a structured categorization of gene activities, which we utilized to examine the shared participation of transcription factors associated with reference or SNP alleles in biological processes directly linked to COVID-19 pathogenesis.33 This tool was used to analyze the overrepresented biological processes characterizing TFs binding to the reference/alternative alleles.
The methods used are summarized in Figure 2.
Results
The analysis identified associations between polymorphic variants rs1043618 HSPA1A, rs13161158 HSPA4, rs910652 HSPA12B, rs7189628 DNAJA2, and rs4926222 DNAJB1 and T2DM in patients with a severe course of COVID-19. rs1043618 HSPA1A (odds ratio [OR] = 1.82, 95% confidence interval [CI] = 1.01–3.27, Pperm = 0.04), rs910652 HSPA12B (OR = 2.23, 95% CI = 1.18–4.21, Pperm = 0.01), and rs4926222 DNAJB1 (OR = 2.15, 95% CI = 1.08–4.27, Pperm = 0.03) were associated with higher odds of T2DM in the total sample, whereas rs13161158 HSPA4 (OR = 0.11, 95% CI = 0.01–0.91, Pperm = 0.0075) and rs7189628 DNAJA2 (OR = 0.08, 95% CI = 0.01–0.65, Pperm = 0.001) were associated with lower odds of T2DM (Table 1).
| Genetic variant | Effect allele | Other allele | N1 | OR [95% CI]2 | Pperm3 |
|---|
| T2DM |
| rs7155973 HSP90AA1 | A | G | 179 | 0.60 [0.16-2.26] | 0.43 |
| rs1043618 HSPA1A | C | G | 185 | 1.82 [1.01-3.27] | 0.04 |
| rs13161158 HSPA4 | C | T | 186 | 0.11 [0.01-0.91] | 0.0075 |
| rs753856 HSPA6 | G | C | 184 | 0.77 [0.32-1.90] | 0.57 |
| rs1461496 HSPA8 | A | G | 186 | 0.68 [0.38-1.24] | 0.2 |
| rs1136141 HSPA8 | A | G | 186 | 1.14 [0.55-2.34] | 0.73 |
| rs10892958 HSPA8 | G | C | 181 | 1.08 [0.53-2.19] | 0.84 |
| rs1042665 HSPA9 | C | T | 185 | 0.93 [0.52-1.69] | 0.82 |
| rs910652 HSPA12B | C | T | 187 | 2.23 [1.18-4.21] | 0.01 |
| rs862832 HSPA12B | T | C | 176 | 1.73 [0.80-3.77] | 0.18 |
| rs17155992 HSPA14 | A | G | 184 | 1.06 [0.39-2.84] | 0.91 |
| rs2034598 DNAJA2 | G | A | 188 | 1.24 [0.68-2.27] | 0.48 |
| rs7189628 DNAJA2 | T | C | 186 | 0.08 [0.01-0.65] | 0.001 |
| rs4926222 DNAJB1 | G | A | 186 | 2.15 [1.08-4.27] | 0.03 |
| rs4279640 HSF1 | C | T | 186 | 0.94 [0.55-1.62] | 0.83 |
| rs6909985 HSF2 | T | G | 184 | 1.14 [0.53-2.43] | 0.74 |
| rs706121 BAG1 | C | T | 185 | 1.56 [0.79-3.11] | 0.21 |
| rs196336 BAG3 | T | C | 185 | 0.82 [0.45-1.50] | 0.53 |
| rs196329 BAG3 | A | G | 182 | 0.77 [0.40-1.50] | 0.44 |
| Obesity |
| rs7155973 HSP90AA1 | A | G | 175 | 0.73 (0.30-1.77) | 0.48 |
| rs1043618 HSPA1A | C | G | 181 | 1.15 (0.74-1.78) | 0.53 |
| rs13161158 HSPA4 | C | T | 182 | 0.74 (0.32-1.72) | 0.48 |
| rs753856 HSPA6 | G | C | 180 | 2.11 (1.07-4.19) | 0.03 |
| rs1461496 HSPA8 | A | G | 182 | 0.94 (0.61-1.45) | 0.78 |
| rs1136141 HSPA8 | A | G | 182 | 1.05 (0.61-1.80) | 0.86 |
| rs10892958 HSPA8 | G | C | 177 | 1.10 (0.65-1.87) | 0.73 |
| rs1042665 HSPA9 | C | T | 181 | 1.30 (0.83-2.04) | 0.24 |
| rs910652 HSPA12B | C | T | 183 | 0.83 (0.52-1.33) | 0.44 |
| rs862832 HSPA12B | T | C | 172 | 1.11 (0.59-2.08) | 0.74 |
| rs17155992 HSPA14 | A | G | 180 | 1.21 (0.57-2.56) | 0.61 |
| rs2034598 DNAJA2 | G | A | 184 | 1.20 (0.76-1.89) | 0.43 |
| rs7189628 DNAJA2 | T | C | 182 | 1.27 (0.62-2.61) | 0.51 |
| rs4926222 DNAJB1 | G | A | 182 | 1.00 (0.57-1.73) | 0.99 |
| rs4279640 HSF1 | C | T | 182 | 0.87 (0.59-1.29) | 0.49 |
| rs6909985 HSF2 | T | G | 180 | 1.06 (0.59-1.90) | 0.84 |
| rs706121 BAG1 | C | T | 181 | 1.21 (0.70-2.10) | 0.49 |
| rs196336 BAG3 | T | C | 181 | 1.33 (0.87-2.05) | 0.19 |
| rs196329 BAG3 | A | G | 178 | 1.06 (0.66-1.72) | 0.81 |
We also observed that rs753856 HSPA6 (effect allele G) was associated with higher odds of obesity in COVID-19 patients (OR = 2.11, 95% CI = 1.07–4.19, Pperm = 0.03) (Table 1). Because no across-SNP multiple-testing correction was applied, all genetic association findings should be considered exploratory pointwise associations requiring independent replication. Because the T2DM subgroup included only 39 patients, the T2DM-related genetic association findings should be interpreted with particular caution.
In exploratory analyses, we observed several nominal associations between the studied SNPs and clinical parameters of COVID-19 patients. We applied a Bonferroni correction (P value multiplied by three) to these results. Nonetheless, given the small genotype-stratified cell counts in several subgroup analyses, these clinical-parameter associations should be regarded as exploratory signals rather than confirmatory findings. In the group of patients without obesity, rs910652 HSPA12B (effect allele C) showed a nominal overall genotype-group difference in fibrinogen levels (Pk-w = 0.04). The C/G genotype of rs753856 HSPA6 in the group with obesity was nominally associated with lower C-reactive protein (CRP) levels (Pbonfm-w = 0.003). Moreover, the rs753856 HSPA6 C/G heterozygotes had more days of oxygen support than C/C homozygotes (Pbonfm-w = 0.0015). The same genotype was nominally associated with lower CRP levels in the group without T2DM (Pbonfm-w = 0.045). SNP rs13161158 HSPA4 (effect allele C) also showed a nominal overall genotype-group difference in fibrinogen levels in the group without T2DM (Pk-w = 0.04) (Supplementary Table. 4).
Functional annotation of cis-eQTL patterns across multiple tissues identified several potentially relevant associations.
Allele C rs1043618 HSPA1A was associated with lower expression of 16 genes: C4A and LY6G5B in subcutaneous adipose, liver, skeletal muscles, whole blood, CYP21A1P and LY6G5C in visceral and subcutaneous adipose tissues, skeletal muscles, and whole blood, DDAH2 in visceral and subcutaneous adipose tissues and skeletal muscles, HLA-DRB5 in visceral and subcutaneous adipose tissues, liver, skeletal muscles, pancreas, and whole blood, ATF6B in visceral adipose and skeletal muscles, NOTCH4 in subcutaneous adipose and skeletal muscles, PRRC2A in skeletal muscles and whole blood, VARS2 in subcutaneous adipose and whole blood, C6orf48, HCG22 in subcutaneous adipose, PBX2 and VWA7 in whole blood, RNF5 and APOM in skeletal muscles.
Moreover, allele C rs1043618 HSPA1A was associated with higher expression of 9 genes: C4B in visceral adipose and whole blood, CCHCR1 in visceral and subcutaneous adipose tissues, CYP21A2 in subcutaneous adipose tissues and whole blood, HLA-DQB2 and MICB in visceral and subcutaneous adipose tissues and skeletal muscles, STK19B in subcutaneous adipose, HSPA1L and HLA-B in pancreas, HLA-C in skeletal muscles.
Allele C of rs13161158 HSPA4 was associated with higher AFF4 expression in subcutaneous adipose tissue and pancreatic tissue.
Allele G rs753856 HSPA6 was associated with higher RPS23P10 expression in whole blood.
The C allele of rs910652 HSPA12B was associated with lower CDC25B expression and higher LINC01730 expression in subcutaneous adipocytes, as well as higher HSPA12B and LINC01730 expression in whole blood (Supplementary Table. 5).
The cis-eQTL profiles of rs1043618 HSPA1A, rs13161158 HSPA4, rs910652 HSPA12B, rs7189628 DNAJA2, and rs4926222 DNAJB1 in peripheral blood are presented in Supplementary Table. 6.
Analysis of transcription-factor (TF) binding predicted several significant differences. The full list of TFs is presented in Supplementary Table. 7, 8, 9, 10, 11.
The risk allele С rs910652 HSPA12B was predicted to create DNA-binding sites for 21 TF, annotated to the following biological processes: retinoic acid receptor signaling pathway (GO:0048384, FDR = 2.86×10-2), cellular response to hypoxia (GO:0071456, FDR = 1.58×10-5), positive regulation of cysteine-type endopeptidase activity involved in apoptotic process (GO:0008656, FDR = 3.22×10-4), glucose homeostasis (GO:0042593, FDR = 9.15×10-5). The protective allele T rs910652 HSPA12B was predicted to create DNA-binding sites for 39 TF, annotated to interleukin-9-mediated signaling pathway (GO:0038113, FDR = 5.57×10-3), interleukin-2-mediated signaling pathway (GO:0038110, FDR = 6.88×10-3), growth hormone receptor signaling pathway via JAK-STAT (GO:0060397, FDR = 8.94×10-3), interleukin-15-mediated signaling pathway (GO:0035723, FDR = 9.9×10-3), cellular response to estrogen stimulus (GO:0071391, FDR = 2.44×10-2).
The risk allele C rs7189628 DNAJA2 was predicted to create DNA-binding sites for 64 TFs annotated to negative regulation of interleukin-4 production (GO:0032713; FDR = 3.69×10-3), negative regulation of interleukin-5 production (GO:0032714; FDR = 4.55×10-3), negative regulation of oxidative stress-induced neuron intrinsic apoptotic signaling pathway (GO:1903377; FDR = 8.00×10-3), positive regulation of vascular endothelial growth factor production (GO:0010575; FDR = 9.05×10-4), regulation of hormone biosynthetic process (GO:0046885; FDR = 2.42×10-2), cytokine production (GO:0001816; FDR = 2.79×10-2), negative regulation of extrinsic apoptotic signaling pathway via death domain receptors (GO:1902042; FDR = 3.04×10-2), B-cell homeostasis (GO:0001782; FDR = 3.19×10-2), regulation of transforming growth factor beta production (GO:0071634; FDR = 3.32×10-2), cellular response to hypoxia (GO:0071456; FDR = 2.47×10-2), transforming growth factor beta receptor superfamily signaling pathway (GO:0141091; FDR = 4.48×10-2), response to oxidative stress (GO:0006979; FDR = 3.33×10-2), cellular response to cytokine stimulus (GO:0071345; FDR = 4.44×10-2). The protective allele T rs7189628 DNAJA2 was predicted to create DNA-binding sites for 28 TFs annotated to response to retinoic acid (GO:0032526; FDR = 1.67×10-2) and cellular response to reactive oxygen species (GO:0034614; FDR = 3.24×10-2).
The risk allele G rs4926222 DNAJB1 was predicted to create DNA-binding sites for 85 TFs annotated to the following biological processes: positive regulation of transcription from RNA polymerase II promoter by glucose (GO:0000432; FDR = 5.11×10-5), peroxisome proliferator-activated receptor signaling pathway (GO:0035357; FDR = 1.35×10-2), response to fructose (GO:0009750; FDR = 1.61×10-2), retinoic acid receptor signaling pathway (GO:0048384; FDR = 3.74×10-2), positive regulation of transforming growth factor beta production (GO:0071636; FDR = 4.10×10-2), cellular response to glucocorticoid stimulus (GO:0071385; FDR = 1.22×10-3), negative regulation of fat cell differentiation (GO:0045599; FDR = 1.85×10-2), steroid hormone receptor signaling pathway (GO:0043401; FDR = 3.71×10-3), cellular response to glucose stimulus (GO:0071333; FDR = 4.46×10-2), response to retinoic acid (GO:0032526; FDR = 1.04×10-2), fat cell differentiation (GO:0045444; FDR = 1.59×10-2), lipid homeostasis (GO:0055088; FDR = 4.82×10-3), regulation of transforming growth factor beta receptor signaling pathway (GO:0017015; FDR = 2.33×10-2), response to hypoxia (GO:0001666; FDR = 8.91×10-4), regulation of hormone secretion (GO:0046883; FDR = 4.04×10-3), response to ketone (GO:1901654; FDR = 1.37×10-2).
The protective allele C rs13161158 HSPA4 was predicted to create DNA-binding sites for 28 TFs annotated to retinoic acid receptor signaling pathway (GO:0048384; FDR = 4.28×10-2). The risk allele T rs13161158 HSPA4 was predicted to create DNA-binding sites for 38 TFs annotated to negative regulation of fat cell proliferation (GO:0070345; FDR = 2.94×10-3), intrinsic apoptotic signaling pathway by p53 class mediator (GO:0072332; FDR = 6.59×10-3), endocrine system development (GO:0035270; FDR = 4.54×10-2), regulation of carbohydrate metabolic process (GO:0006109; FDR = 7.46×10-3), cellular response to oxygen-containing compound (GO:1901701; FDR = 1.64×10-3), response to hormone (GO:0009725; FDR = 2.22×10-2).
The protective allele C rs753856 HSPA6 was predicted to create DNA-binding sites for 39 TFs annotated to the following biological processes: interleukin-4-mediated signaling pathway (GO:0035771; FDR = 5.28×10-3), growth hormone receptor signaling pathway via JAK-STAT (GO:0060397; FDR = 9.52×10-3), cellular response to interleukin-6 (GO:0071354; FDR = 1.16×10-3), non-canonical NF-kappaB signal transduction (GO:0038061; FDR = 1.51×10-3), canonical NF-kappaB signal transduction (GO:0007249; FDR = 4.12×10-3), cellular response to hydrogen peroxide (GO:0070301; FDR = 1.01×10-2) (Fig. 3).
Discussion
To our knowledge, this study provides preliminary evidence that rs1043618 HSPA1A, rs910652 HSPA12B, and rs4926222 DNAJB1 were associated with higher odds of T2DM among patients with severe COVID-19, whereas rs13161158 HSPA4 and rs7189628 DNAJA2 were associated with lower odds of T2DM. Interestingly, rs1043618 has previously been associated with T2DM, which is consistent with our results.34 Furthermore, rs753856 HSPA6 was associated with higher odds of obesity in patients with severe COVID-19.
HSPA1A, HSPA12B, HSPA4, and HSPA6 encode members of the 70-kDa HSP family (HSP70s), which are involved in the protein quality control system.35 These chaperones participate in the ubiquitin-proteasome system, which is one of the major pathways regulating proteostasis.36 Our cis-eQTL analysis suggests that SNP rs1043618 HSPA1A may influence expression of the RNF5 gene, which encodes an E3 ubiquitin-protein ligase.
Furthermore, HSP70s can interrupt nuclear factor kappa B (NF-κB) signaling by inhibiting ubiquitination of key signaling transducers.23 NF-κB signaling plays an essential role in inflammatory response via inducing synthesis of proinflammatory cytokines.37 Recent studies have associated NF-κB hyperactivation in patients with COVID-19 with dysregulated immune responses and cytokine storm.38,39 Notably, NF-κB signaling is enhanced in patients with metabolic syndrome, making them more susceptible to critical illness due to COVID-19.40,41 NF-κB signaling is also involved in low-grade obesity-associated inflammation: overnutrition triggers hypothalamic IKKβ/NF-κB activation, which disrupts insulin/leptin signaling.42,43 Indeed, we have observed that the sequence surrounding rs753856 HSPA6 was predicted to bind several TFs involved in NF-κB signal transduction.
Endoplasmic reticulum (ER) stress has been implicated as an inducer of hypothalamic IKKβ/NF-kappaB pathway.43 Prolonged ER stress responses have been observed in cells exposed to hyperglycemia and hyperlipidemia and may contribute to cytokine-induced β-cell death.44 Elevated ER stress in pancreatic cells has been associated with cellular dysfunction and subsequent insulin deficiency.45,46 An association between increased ER stress and HSP70 overexpression has been reported in several cancer cell models.47,48 Notably, we observed that rs1043618 HSPA1A may be related to ER stress pathways through cis-eQTL-mediated differences in ATF6B expression. ATF6B encodes a transcription factor that activates unfolded protein response target genes during ER stress.49
Members of the 40-kDa HSP family encoded by DNAJA2 and DNAJB1 genes, are also involved in unfolded protein response. Several members of the DnaJ protein family may mitigate ER stress in hyperglycaemic conditions.50,51DNAJA2 provides cellular protection against reactive oxygen species.52 Elevated glucose levels increase production of reactive oxygen species, subsequently promoting pancreatic β-cell apoptosis.53 Therefore, effective defense against oxidative stress may help protect β-cells from oxidative injury.54 Our TF analysis suggests that the T allele of rs7189628 DNAJA2 is involved in cellular response to reactive oxygen species.
DNAJB1 serves as HSPA1A co-chaperone, assisting in proteostasis maintenance.55 Levels of DNAJB1 in patients with T2DM are consistently higher than in controls, suggesting that its expression is upregulated under metabolic stress.55,56HSPA1A itself participates in antigen presentation, stimulates innate and adaptive immune responses,19 regulates cytokine production.57 Cis-eQTL analysis suggests that rs1043618 HSPA1A may influence expression of genes encoding human leukocyte antigens (HLA): HLA-DRB5, HLA-DQB2, HLA-C, and HLA-B. Expression of the MICB gene may also be influenced by rs1043618 HSPA1A. MICB encodes a protein involved in major histocompatibility complex class I antigen processing and presentation.58 It promotes cytotoxicity and cytokine production under conditions of viral infection.59
It is important to note that rs1043618 HSPA1A may influence expression of the PRRC2A and DDAH2 genes, which are linked to T2DM pathogenesis. Several genomic studies have reported associations between PRRC2A and obesity, insulin resistance, and T2DM.60,61 Moreover, A PRRC2A-dependent mechanism has been reported to regulate inflammation and glycolysis in macrophages via m6A/Npas2/HIF-1α axis in experimental models of diabetic nephropathy,62 thus participating in T2DM-associated inflammation. Increased asymmetric dimethylarginine (ADMA) levels are characteristic of T2DM patients.63DDAH2 has been reported to enhance glucose-dependent insulin secretion and DDAH1/DDAH2 genetic variation is associated with ADMA levels in diabetic individuals.64,65 ADMA concentrations are associated with T2DM complications.63,66 Furthermore, high DDAH2 levels have been reported to suppress antiviral immune responses.67
Among the studied polymorphisms in genes encoding HSP70s, only rs13161158 HSPA4 was associated with lower odds of T2DM. Our in silico analysis predicted that this locus alters DNA-binding sites for TFs annotated to negative regulation of fat cell proliferation. These results should be interpreted cautiously. Supplementary Table. 5 and 6 indicate that rs13161158 HSPA4 is associated with increased AFF4 expression; because AFF4 is a mediator of adipogenic differentiation, this finding may provide a possible link between this locus and metabolic regulation.68
Limitations
Several limitations should be acknowledged. First, this investigation was conducted in a two-center cohort in a Central Russian population and included exclusively ICU patients. It lacked an independent replication cohort; therefore, the generalizability of the study to other ethnic groups and to patients with mild or moderate COVID-19 or non-ICU populations is limited. In addition, genetic ancestry or population structure was not formally evaluated; therefore, residual population stratification cannot be excluded. Second, the absence of a formal pre-study power calculation and the small T2DM subgroup limits statistical power and may make effect estimates unstable. Third, sparse genotype-stratified analyses increase uncertainty and residual confounding remains possible despite covariate adjustment. Fourth, gene expression and protein levels were not measured in HSP genotype carriers; hence, the proposed molecular mechanisms rely on in silico functional annotation and require experimental validation. Finally, the exploratory analyses of clinical parameters involved multiple SNPs, subgroups, and clinical outcomes, and should be interpreted cautiously until replicated in larger independent cohorts. All findings are associations rather than causal effects.
Conclusions
This study analyzed associations between HSP SNPs and the risk of T2DM and obesity among patients with severe COVID-19. It was an exploratory case-control study among ICU patients with severe or critical COVID-19 recruited in Central Russia. We also observed exploratory nominal associations between HSP SNPs and clinical parameters such as thrombodynamic measures, duration of oxygen therapy, and CRP levels. Functional annotation suggests plausible mechanisms by which risk-associated HSP SNPs may be involved in T2DM and obesity, but these findings require independent replication and functional validation. The findings represent associations rather than causal effects and should be interpreted in light of the small T2DM subgroup and the absence of a pre-study power calculation.
Supporting information
Supplementary material for this article is available at https://doi.org/10.14218/GE.2026.00015 .
Supplementary Table. 1
Baseline and clinical characteristics stratified by T2DM status
(DOCX)
Supplementary Table. 2
Baseline and clinical characteristics stratified by obesity status
(DOCX)
Supplementary Table. 3
Quality-control metrics and SNP-specific regression sample sizes for the genotyped HSP SNPs
(DOCX)
Supplementary Table. 4
Clinical and laboratory parameters of COVID-19 patients (with/without T2DM or obesity)
(DOCX)
Supplementary Table. 7
Analysis of the effect of rs910652 HSPA12B on the binding of DNA to transcription factors
(DOCX)
Supplementary Table. 8
Analysis of the effect of rs7189628 DNAJA2 on the binding of DNA to transcription factors
(DOCX)
Supplementary Table. 9
Analysis of the effect of rs4926222 DNAJB1 on the binding of DNA to transcription factors
(DOCX)
Supplementary Table. 10
Analysis of the effect of rs13161158 HSPA4 on the binding of DNA to transcription factors
(DOCX)
Supplementary Table. 11
Analysis of the effect of rs753856 HSPA6 on the binding of DNA to transcription factors
(DOCX)
Declarations
Funding
The study was funded by Kursk State Medical University.
Conflict of interest
The authors declare no conflicts of interest related to this publication.
Author contributions
Study concept and design (KK, OB), acquisition of data (AK, VS), analysis and interpretation of data (KK, EL, AA), drafting of the manuscript (KK, AD), critical revision of the manuscript for important intellectual content (OB), administrative, technical, or material support (OB), and study supervision (OB). All authors have made significant contributions to this study and have approved the final manuscript.
Ethical statement
The study was conducted in accordance with the Declaration of Helsinki (as revised in 2024), and was approved by the Ethical Review Committee of Kursk State Medical University (protocol No. 1, January 11, 2022). Written informed consent was obtained from all individual participants included in the study.
Data sharing statement
Individual-level genotype and clinical data are not publicly available because of ethical and privacy restrictions. Aggregated data supporting the main findings are provided in the manuscript and supplementary materials.