Introduction
Metabolic syndrome (MetS) is a cluster of metabolic conditions characterized by obesity, insulin resistance (IR), hyperglycemia, hypertension, and dyslipidemia.1 Accumulated data now indicate that MetS is considered a chronic inflammatory disease.2 One of the most important molecular links between immunologic and metabolic processes in MetS is mitochondrial dysfunction, characterized by disruption of mitochondrial spatial integrity and genome, increased production of reactive oxygen species, and decreased ATP production.3 Mitochondrial dysfunction in immune cells mediates meta-inflammation.4 Thus, the state of mitochondria can determine the transition from an anti-inflammatory to a pro-inflammatory phenotype and regulate cytokine secretion.5 Specific manifestations of mitochondrial dysfunction in MetS include decreased copy number and accumulation of mutations in mitochondrial DNA (mtDNA).6 In addition, extracellular mtDNA acts as a damage-associated molecular pattern and contributes to the development and maintenance of chronic inflammation.3,5
It is important to note that mtDNA is a multicopy molecule—unlike nuclear DNA, which contains one (haploid) or two (diploid) copies of the genome per cell, mtDNA is present in multiple copies per cell, ranging from hundreds to thousands depending on the tissue type.7 Not all of these copies are identical; therefore, mtDNA mutations can have varying degrees of genetic expressivity and contribute differently to the development of mitochondrial dysfunction. This condition is called mitochondrial heteroplasmy.8 It is believed that up to a certain “phenotypic threshold” of heteroplasmy, the cell retains the ability to compensate for mutant mtDNA.9 However, the “phenotypic threshold” is highly specific to the cell type and its state and also depends on the nature of the genetic variant (pathological, neutral, or positive).9 Therefore, determining the “normal” level of heteroplasmy is very difficult.
In general, evolution opposes mitochondrial heteroplasmy—the norm for eukaryotes is maternal gene transmission, and male mtDNA is actively eliminated during fertilization. A longitudinal study by Lechuga-Vieco AV et al. (2022) on the creation of heteroplasmic mice (fusion of two normal mitochondrial haplotypes at the zygote stage) showed that even non-pathological heteroplasmy causes serious consequences for tissues unable to resolve heteroplasmy, such as cardiac and skeletal muscles and the eyes. Non-pathological heteroplasmy-induced metabolic stress in these animals leads to severe pathology in adulthood, including pulmonary hypertension, severe heart failure, skeletal muscle wasting, weakness, and an increased risk of premature death.10
While the effect of heteroplasmy on the phenotype of cells in muscle, cardiac, and vascular tissues is being actively studied, the effect of heteroplasmy on the metabolism of immune cells, which play a significant role in the development of MetS (as discussed above), remains unexplored. Publications indicate that blood monocytes are particularly promising for study in the context of MetS pathogenesis because:
They are among the most highly heteroplasmic populations of immune cells in the blood. A study by Walker et al.11 showed that three patients with MELAS (mitochondrial encephalomyopathy with lactic acidosis and stroke-like episodes) had the highest levels of heteroplasmy at position m.3243A>G (the most common pathogenic heteroplasmic mtDNA mutation), detected specifically in monocytes and dendritic cells. In another study, mtDNA from monocytes of healthy donors, compared with mtDNA from B lymphocytes, had a higher number of deleterious nonsynonymous substitutions (according to CADD score).12
Monocytes actively participate in the pathogenesis of MetS. In a “bad” cardiometabolic environment, they acquire a pro-inflammatory phenotype and are recruited to adipose tissue and the liver, where they differentiate into macrophages and contribute to the formation of a pro-inflammatory micro- and macroenvironment.13
Several studies report a significant role of heteroplasmy in monocytes in various MetS components,14–16 such as obesity and cardiovascular diseases. However, most of these studies analyze only specific heteroplasmic points; they do not screen for genetic, molecular, or biochemical analytes and do not attempt to establish an association between functional indicators and detectable heteroplasmy.
Therefore, this study aimed to identify associations of monocytic mtDNA variability with its phenotypic indices, including cytokine secretion and gene expression, and cardiometabolic parameters in patients with MetS.
Materials and methods
Study design
The study design is schematically presented in Figure 1.
Study participants
The diagnosis of the study participants was verified by an endocrinologist at the Clinical Diagnostic Center of the Immanuel Kant Baltic Federal University. The study was conducted in accordance with the World Medical Association Declaration of Helsinki (2024) and the Protocol to the Convention on Human Rights and Biomedicine (1999). The ethics committee of Immanuel Kant Baltic Federal University approved the study (IKBFU Ethics Committee Conclusion No. 44, dated 13.12.2023). All participants were selected and recruited for the study by an endocrinologist, and each provided signed informed consent.
Venous blood was drawn from the patients by puncturing the cubital vein and was collected in the morning while fasting. Each patient was assigned an internal serial number for anonymous data processing.
Blood biochemistry tests and cardiometabolic parameters
Blood biochemistry tests and triplex scanning were performed to assess the patients’ cardiometabolic profiles. Blood samples for biochemical analysis were collected in 5 mL tubes containing silicon dioxide as a clot activator. Plasma biochemistry was analyzed using a Furuno CA-180 biochemistry analyzer (Furuno Electric Company, Nishinomiya City, Japan) and DiaSys test systems (DiaSys Diagnostic Systems, Holzheim, Germany). Insulin levels were determined using an enzyme-linked immunosorbent assay (ELISA) (Vector-Best, Novosibirsk, Russia). The HOMA-IR (Homeostatic Model Assessment of Insulin Resistance) index was calculated using the formula:
Insulin (mIU/L) × Fasting blood glucose (mmol/L)22.5.
The atherogenic coefficient was calculated using the formula:
Total cholesterol (mmol/L)-High-density lipoprotein (mmol/L)High-density lipoprotein (mmol/L).
Triplex scanning of the neck veins was performed by an ultrasound physician at the recommendation of an endocrinologist. Intima-media thickness and percentage of stenosis were determined. Descriptive statistics for the blood biochemistry test are presented as median (Q1; Q3). Comparisons between groups were made using one-way analysis of variance (ANOVA) with the Tukey comparison test or the Kruskal-Wallis test with Dunn’s multiple comparisons, depending on the data distribution, which was assessed by the Shapiro-Wilk test.
Monocyte separation and culturing conditions
Monocytes were isolated to analyze cytokine secretion in vitro. Blood for CD14+ cell separation was collected in a volume of 30 mL using vacutainer tubes containing EDTA anticoagulant. Peripheral blood mononuclear cells were isolated from whole blood using a Ficoll gradient (1.077 g/mL) (Dia-Em, Moscow, Russia). CD14+ cells were then isolated from the peripheral blood mononuclear cells by immunomagnetic separation. A total of 100 µL of CD14+ magnetic beads (Miltenyi Biotec, Bergisch Gladbach, Germany) were added to the cell suspension. Separation was performed using a MACS Multi system and LS Separation Columns. CD14+ cells were counted with a Countess Automated Cell Counter (Invitrogen, Waltham, Massachusetts, USA) using 0.4% Trypan Blue (Invitrogen, Waltham, Massachusetts, USA).
After several washes with PBS, the desired number of separated CD14+ cells was resuspended in Iscove’s serum-free growth medium (Gibco, Waltham, Massachusetts, USA) and added at a rate of 5 × 105 cells to two wells of a 48-well plate for each donor. Lipopolysaccharide (LPS) (Sigma, St. Louis, Missouri, USA) was added to one well to a final concentration of 0.1 µg/mL. The cells were cultivated in 5% CO2 at 37 °C for 24 h. After 24 h of culturing, the growth medium was changed, and the supernatant was collected and stored at −80 °C. A total of 1 × 106 separated CD14+ cells in Iscove’s medium was collected for subsequent DNA extraction and stored at −20 °C.
All remaining cells from the initial population of separated CD14+ cells were centrifuged at 8,000 g for 5 min. After centrifugation, the supernatant was removed, and the pellet was resuspended in 350 µL of ExtractRNA (Evrogen, Moscow, Russia), incubated for 10 min at room temperature, and then stored at −80 °C until subsequent RNA extraction.
To assess cytokine secretion (ELISA), expression (quantitative real-time polymerase chain reaction (RT-qPCR)), and variability of mtDNA (next-generation sequencing) in human monocytes (CD14+ cells), the following samples were obtained from each patient:
1 × 106 cells in Iscove’s medium for DNA extraction – native CD14+ cells;
a variable number of cells in ExtractRNA for RNA extraction – native CD14+ cells;
cell supernatant after 24 h of culture without LPS for cytokine level measurement – basal secretion of CD14+ cells;
cell supernatant after 24 h of culture with LPS for cytokine level measurement – stimulated secretion of CD14+ cells.
mtDNA sequencing
For the mtDNA variability study, mtDNA sequencing was performed using next-generation sequencing technology for 87 subjects in our study. Total DNA was isolated from human CD14+ cells using the ExtractDNA Blood & Cells total DNA isolation kit (Evrogen, Moscow, Russia) according to the manufacturer’s protocol.
Full-length mtDNA sequencing was conducted in several stages. Long-range polymerase chain reaction (PCR) was conducted using the BioMaster LR HS-PCR (2×) amplification kit (BioLabMix, Novosibirsk, Russia). PCR was run on a 2720 Thermal Cycler (Applied Biosystems, Foster City, California, USA) according to the manufacturer’s protocol. The quality of the PCR products was assessed by electrophoresis in a 1.5% agarose gel. The concentration of the PCR products was measured on a Qubit fluorometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA) using the Spectra Q BR DNA concentration reagent kit (Raissol, Moscow, Russia). PCR product mixtures for each sample were used to prepare DNA libraries with the SG GM Plus whole-genome library reagent kit (Raissol, Moscow, Russia). Sample preparation was performed according to the manufacturer’s protocol without modification. The quality of the libraries was checked on a TapeStation 4200 instrument (Agilent, Santa Clara, California, USA). The concentration of each sample was then determined on a Qubit fluorometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA) using the Spectra Q BR or Spectra Q HS reagent kit (Raissol, Moscow, Russia). DNA library sequencing was performed on a GenoLab M sequencer (GeneMind, Luohu, Shenzhen, China) using the GenoLab M V2.0 sequencing kit (FCM cell, 150 cycles, 250M reads per cell) (GeneMind, Luohu, Shenzhen, China). Mean sequencing coverage was 10,000, and reads were paired-end with a length of 75 bp. FASTQ files were aligned to the human genome assembly hg38 using DRAGEN 3.9.5 DNA pipeline software (Illumina, San Diego, California, USA).
To identify positions with homoplasmy and heteroplasmy, the number of reads for each of the four nucleotides at each genome position related to chrM was determined. Sequencing artifacts were analyzed as follows:
when strand bias was detected (a significant skew in the number of reads from one DNA strand compared to the other), heteroplasmy was not considered reliable;
heteroplasmies located in the primer annealing zone were also not considered reliable;
heteroplasmy below 5% or above 95% was not considered reliable.
Significant mutations in mtDNA were classified as follows using the approach described in Calabrese et al.15: low heteroplasmy (positions where 5% to 10% of reads were alternative nucleotides relative to the reference genome), intermediate heteroplasmy (10–95% heteroplasmy), and alternative homoplasmy (≥95%). Calabrese et al. showed that using these thresholds (5–10% and 10–95%) allows for the delineation of heteroplasmies with more prominent pathological effects. mtDNA variability was quantitatively analyzed using the parameters of level, number, prevalence, and frequency of mitochondrial single nucleotide variants (mtSNVs). Definitions for these terms are provided in Table 1. The data on the amount of mtSNVs per individual were presented as mean ± standard deviation. P-values were calculated using one-way ANOVA with the Tukey comparison test or the Kruskal-Wallis test with Dunn’s multiple comparisons, depending on the data distribution, which was assessed by the Shapiro-Wilk test.
Table 1Terms and definitions of mtDNA variability parameters
| Term | Definition | Calculation and units of measurement |
|---|
| MtSNV level | The proportion of reads containing the alternative nucleotide at the position | The number of reads of alternative nucleotides at a position divided by the total read depth at that position – Proportion or percentage |
| Amount of mtSNVs per individual | The number of reads of alternative nucleotides at a position divided by the total read depth at that position | Counting – Absolute value |
| MtSNV prevalence | The proportion of people in a group who have this mtSNV | The number of people with the particular mtSNV divided by the total number of people in the group – Proportion or percentage |
| MtSNV frequency | Frequency of occurrence of mtSNV at a specific locus (entire mtDNA or an individual locus) | Number of mtSNVs at a locus divided by the length of the locus – Proportion or percentage |
Gene expression analysis
For gene expression analysis, RNA from human monocytes was extracted, and RT-qPCR was performed. Total RNA was isolated using the ExtractRNA reagent (Evrogen, Moscow, Russia) based on the phenol–chloroform method, according to the manufacturer’s protocol. The total RNA was eluted in 50 µL of RNase-free water. RNA concentration was measured immediately after extraction using an Implen NanoPhotometer N (Implen, München, Germany). Samples were stored at −80 °C until reverse transcription, followed by RT-qPCR.
RT-qPCR for measuring gene expression levels was performed as follows: universal reverse transcription to obtain cDNA from mRNA of human CD14+ cells was carried out using the MMLV RT kit (Evrogen, Moscow, Russia) on equimolar volumes of RNA, with RiboCare RNase inhibitor (Evrogen, Moscow, Russia) added to increase reaction efficiency. Otherwise, reverse transcription followed the manufacturer’s protocol. For the negative reverse transcription control, an equal volume of RNase-free water was added instead of RNA.
The cDNA generated from the reverse transcription reaction was diluted 26-fold and used for quantitative PCR with the 5× qPCRmix-HS SYBR PCR mixture (Evrogen, Moscow, Russia), according to the manufacturer’s protocol. Primers were added to a final concentration of 0.4 µM (synthesized by BioBeagle, Saint Petersburg, Russia). Four groups of genes were studied: 1) mitochondrial biogenesis genes—mtDNA transcription factors (TFAM, NFE2L2) and genes responsible for mitochondrial division and fusion (MFN2, OPA1, DRP1); 2) mitochondrial uncoupling genes, a process that reduces the membrane potential of mitochondria—UCP2, ANT2; 3) oxidative stress system genes—oxidizer NOX2, antioxidant SOD2; 4) main NF-κB signaling genes—p65, NFKBIA. Primer sequences are provided in Supplementary Table 1. Amplification and real-time PCR data acquisition were performed on a CFX96 thermal cycler (Bio-Rad, Hercules, California, USA), with initial activation at 95 °C for 5 min, followed by 40 cycles of denaturation at 95 °C for 15 s, annealing at the optimal primer temperature for 30 s, and elongation at 72 °C for 30 s. Optimal annealing temperatures were determined empirically using gradient PCR. After amplification, melting curve analysis was performed to confirm reaction specificity.
PCR was performed according to the MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments) guidelines.17 The experiment was designed using sample maximization and included a calibration curve (five standard points in duplicate) for each analyte gene to determine reaction efficiency,17,18 as well as inter-run calibrators.18 For the calibration curve standard points, a pre-generated pool of cDNA collected from all studied samples at a concentration of 0.16 ng/µL was used. The dilution factor was 2.0 (cDNA per reaction at points 1–5: 0.5, 0.25, 0.125, 0.0625, and 0.0313 ng, respectively). The cDNA pool at 0.16 ng/µL was also used as the template for the inter-run calibrators. Expression analysis of the genes of interest was carried out in duplicate for each sample (0.5 ng cDNA per reaction). Each run included a no-template control and a negative control.
All results with different amplicon melting temperatures and a threshold cycle higher than that of the no-template control were excluded from further calculations. To normalize the expression data, a panel of four reference genes (RPL13A, PPIB, PPIA, ACTB) was used; their stability as a reference system was determined using the geNorm protocol. The qBase+ program was used to calculate the PCR results.18
The data were log-transformed to base 10 to normalize the distribution. The transformed expression profile data were tested for normality using the Shapiro–Wilk test. Equality of sample means was then tested using a standard one-way analysis of variance with Tukey’s post hoc test or, for non-normally distributed data, the Kruskal–Wallis test with Dunn’s post hoc test. Differences were considered significant at a P-value less than 0.05. Statistical analysis was performed using GraphPad Prism 9.3.1. Descriptive statistics and group comparisons of gene expression levels are presented in Supplementary Table 2.
Cytokine secretion analysis
ELISA was performed to analyze cytokine secretion, after which the stimulation index was calculated. After 24 h of culturing, the resulting samples (cell culture supernatants) were analyzed for levels of pro- and anti-inflammatory cytokines (interleukin (IL)-6, monocyte chemoattractant protein-1, IL-1β, tumor necrosis factor, IL-10, IL-8) using a three-step sandwich ELISA (all kits from Vector-Best, Novosibirsk, Russia).
Statistical data processing was conducted as follows: the cytokine secretion stimulation index was calculated as the fold change (FC), defined as LPS-stimulated secretion divided by basal secretion without added LPS. FC values were log-transformed to base 2 for normalization. P-values were calculated using one-way ANOVA with the Tukey comparison test or the Kruskal–Wallis test with Dunn’s multiple comparisons, depending on the data distribution, which was assessed by the Shapiro–Wilk test. A P-value less than 0.05 was considered significant. Descriptive statistics for the cytokine secretion stimulation index and the results of group comparisons for these parameters are presented in Supplementary Table 3 (data presented as Median (Q1; Q3)). Statistical analysis was performed using GraphPad Prism 9.3.1.
Statistical analysis of the impact of mtSNV frequency on phenotypic traits
To analyze associations between mtDNA variability and cytokine or gene expression profiles, various statistical tests were performed. Spearman’s rho was used to assess the association between mtSNV frequencies and cardiometabolic parameters. Confounding factors (sex, age, BMI) were excluded by partial Spearman correlation, and associations with a P-value less than 0.05 were considered statistically significant. The data were presented as Median (Q1; Q3).
Multivariate analysis of variance was used to assess the association between mtSNV frequencies, gene expression levels, and cytokine secretion. The primary objective was to determine the extent to which variations in the expression of functional gene groups and cytokine profiles could be explained by mtSNV frequencies in various mitochondrial genome loci (D-loop, protein-coding, tRNA, rRNA, and entire mtDNA).
The studied genes were divided into the following functional groups: 1) mitochondrial biogenesis genes—mtDNA transcription factors (TFAM, NFE2L2) and genes responsible for mitochondrial division and fusion (MFN2, OPA1, DRP1); 2) mitochondrial uncoupling genes, a process that reduces the membrane potential of mitochondria—UCP2, ANT2; 3) oxidative stress system genes—oxidizer NOX2, antioxidant SOD2; 4) main NF-κB signaling genes—p65, NFKBIA.
For each set of dependent variables (cytokines or expression of functional gene groups), a Bray–Curtis pairwise distance matrix was calculated to reflect differences between patients based on parameter profiles. Principal coordinates analysis was performed to represent multivariate data as several orthogonal components that preserve most of the variance. The data distribution type was determined using the Henze–Zirkler or Royston test. Before PERMANOVA (permutational MANOVA), homogeneity of variances between groups was assessed using the betadisper function, followed by PERMANOVA based on 999 permutations. Post hoc pairwise comparisons (pairwise PERMANOVA/MANOVA) were conducted to identify differences between individual groups. A Bonferroni correction for multiple comparisons was applied: a total of 15 comparisons were performed for each set of dependent variables, so only results with a P-value < 0.0033 (0.05/15) were considered significant. The analysis was performed using the Python programming language, version 3.10, with the specified libraries: pandas (2.0) and numpy (1.24) for data preparation and processing, scipy (1.10) for statistical procedures, scikit-bio (0.5) for calculating Bray–Curtis distances, performing PERMANOVA, and assessing homogeneity of variances, pingouin (0.5) for multivariate normality tests (Henze–Zirkler/Royston), and scikit-learn (1.3) for implementing principal coordinate analysis.
Results
Clinical, biochemical and cardiometabolic characteristics
This cross-sectional study recruited three groups of participants, including 1) conditionally healthy blood donors (Control) with a BMI < 25 kg/m2, no abnormalities in cardiometabolic parameters, and no MetS components according to the 2006 International Diabetes Federation (IDF) MetS criteria19; 2) obese patients with a BMI > 30 kg/m2 and no more than one MetS component (according to the 2006 IDF criteria) (Obesity); and 3) patients with a BMI > 30 kg/m2 and at least two MetS components (according to the 2006 IDF criteria) (MetS). A sample of 87 participants had a 75% confidence interval and was calculated by a special formula with a population proportion of 20%.20 Inclusion criteria were: age between 18 and 65 years; absence of any cardiometabolic changes, as determined by an endocrinologist, for the control group; and presence of cardiometabolic changes, as determined by an endocrinologist according to the IDF 2006 criteria. Participants under 18 or over 65, those with acute or severe chronic somatic or infectious diseases, or malignant neoplasms were excluded from the study. The cardiometabolic group characteristics are presented in Table 2.
Table 2Cardiometabolic characteristics of the groups studied
| Variables | Control (1) | Obesity (2) | MetS (3) | P-value | Test |
|---|
| Sex (men/women) | 7 / 27 | 5 / 16 | 11 / 21 | 0.4254 | c |
| Median (Q1; Q3) | 1 vs 2 | 1 vs 3 | 2 vs 3 | |
| Age, years | 39.5 (33.7; 46.0) | 39.0 (33.0; 48.0) | 44.0 (36.0; 49.0) | 0.8241 | 0.2495 | 0.6913 | a |
| BMI, kg/m2 | 23.21 (21.85; 25.19) | 35.70 (31.25; 40.35) | 36.08 (33.66; 41.03) | <0.0001 | <0.0001 | >0.9999 | b |
| Intima media thickness, mm | 0.66 (0.60; 0.72) | 0.70 (0.60; 0.76) | 0.75 (0.60; 0.80) | 0.9423 | 0.7407 | 0.6374 | a |
| Stenosis, % | 0 (0; 0) | 0 (0; 5) | 0 (0; 0) | >0.9999 | >0.9999 | >0.9999 | c |
| Alanine aminotransferase (ALAT), u/L | 12.05 (8.63; 17.23) | 22.40 (16.05; 37.25) | 26.85 (15.05; 41.68) | 0.0001 | <0.0001 | >0.9999 | b |
| Aspartate aminotransferase (ASAT), u/L | 17.10 (15.95; 20.10) | 19.20 (17.90; 25.75) | 22.35 (17.83; 31.70) | 0.033 | 0.002 | >0.9999 | b |
| Fasting blood glucose (FBG), mmol/L | 4.48 (4.32; 4.65) | 4.980 (4.79; 5.31) | 5.21 (4.66; 5.88) | <0.0001 | <0.0001 | >0.9999 | b |
| HOMA-IR | 0.59 (0.55; 0.75) | 1.34 (1.11; 2.23) | 1.70 (1.32; 3.15) | <0.0001 | <0.0001 | 0.8467 | b |
| Insulin, mIU/L | 2.94 (2.77; 3.52) | 5.90 (5.10; 9.31) | 7.26 (5.78; 12.70) | <0.0001 | <0.0001 | 0.7153 | b |
| Triglycerides (TG), mmol/L | 0.69 (0.58; 0.85) | 1.10 (0.84; 1.33) | 1.99 (1.35; 2.37) | 0.0051 | <0.0001 | 0.0276 | b |
| High density lipoprotein (HDL), mmol/L | 1.95 (1.64; 2.05) | 1.47 (1.35; 1.69) | 1.26 (1.01; 1.45) | 0.0001 | <0.0001 | 0.052 | a |
| Low density lipoprotein (LDL), mmol/L | 2.92 (2.60; 3.27) | 2.95 (2.49; 3.34) | 2.86 (2.43; 3.30) | 0.8498 | 0.9718 | 0.9391 | a |
| Non-HDL, mmol/L | 3.70 (2.92; 4.06) | 3.80 (3.07; 4.31) | 4.04 (3.39; 4.59) | 0.5809 | 0.0505 | 0.5144 | a |
| Atherogenicity coefficient | 1.93 (1.58; 2.30) | 2.73 (2.06; 3.09) | 3.19 (2.52; 4.28) | 0.0073 | <0.0001 | 0.0116 | a |
| Total cholesterol (TC), mmol/L | 5.42 (4.71; 5.96) | 5.48 (4.46; 6.07) | 5.36 (4.71; 5.84) | 0.8889 | 0.974 | 0.962 | a |
| P-amylase, u/L | 18.05 (12.43; 23.48) | 11.00 (5.15; 17.55) | 11.90 (7.28; 19.98) | 0.047 | 0.2864 | 0.5349 | a |
| Gamma-glutamyl transferase (GGT), u/L | 12.15 (10.60; 17.70) | 26.20 (18.45; 38.25) | 31.50 (14.23; 62.25) | 0.0001 | <0.0001 | >0.9999 | b |
| Alkaline phosphatase (ALP), u/L | 127.00 (103.00; 161.00) | 185.50 (142.30; 217.30) | 186.50 (154.50; 209.30) | 0.0029 | 0.0007 | >0.9999 | b |
| Uric acid (UA), µmol/L | 263.00 (218.50; 290.20) | 332.90 (281.80; 452.30) | 374.30 (297.00; 444.40) | 0.0005 | <0.0001 | >0.9999 | b |
| Blood urea nitrogen (BUN), mmol/L | 5.30 (4.88; 5.93) | 6.10 (5.15; 6.40) | 5.85 (5.30; 6.70) | 0.2185 | 0.0874 | 0.9741 | a |
| Creatinine (CREA), µmol/L | 72.00 (67.00; 81.25) | 74.00 (69.00; 81.00) | 77.50 (72.00; 84.50) | 0.7984 | 0.134 | 0.5362 | a |
| Total bilirubin (BIL T), µmol/L | 14.48 (12.53; 18.65) | 11.66 (10.28; 13.83) | 11.98 (9.11; 16.54) | 0.0173 | 0.0537 | >0.9999 | b |
| Direct bilirubin (BIL D), µmol/L | 2.79 (2.31; 3.98) | 2.49 (2.03; 2.84) | 2.39 (1.98; 3.32) | 0.1912 | 0.2932 | >0.9999 | b |
| Indirect bilirubin (BIL IND), µmol/L | 11.71 (10.08; 14.37) | 8.88 (8.13; 11.03) | 9.45 (7.12; 12.56) | 0.0132 | 0.0444 | >0.9999 | b |
| Total protein (TP), g/L | 72.80 (70.68; 74.90) | 75.40 (73.55; 76.90) | 75.05 (73.23; 76.88) | 0.0087 | 0.0048 | 0.9852 | a |
| C-reactive protein (CRP), mg/L | 0.26 (0.16; 0.55) | 3.78 (2.05; 5.49) | 5.55 (2.08; 11.40) | <0.0001 | <0.0001 | >0.9999 | b |
| Calcium, mmol/L | 2.37 (2.32; 2.46) | 2.35 (2.27; 2.41) | 2.30 (2.23; 2.46) | 0.7149 | 0.5114 | 0.9628 | a |
| Iron, µmol/L | 16.40 (12.50; 22.40) | 16.80 (12.15; 20.13) | 12.30 (9.10; 18.78) | >0.9999 | 0.2298 | 0.7638 | b |
The amount of heteroplasmic mtSNVs in human CD14+ cells is very low across all study groups
Next-generation sequencing of mtDNA in human CD14+ cells was performed in 87 individuals. The following mtDNA variability parameters were analyzed (Table 2): amount of mtSNVs per individual, mtSNV level, mtSNV prevalence, and mtSNV frequency. mtSNVs were classified into three groups based on mtSNV level: low heteroplasmic (5–10%), intermediate heteroplasmic (10–95%), and alternative homoplasmic mtSNVs (≥95%).
The most common mtSNVs were alternative homoplasmies (all data are presented as mean ± standard deviation of mtSNV amount per individual): Control group 23.11 ± 10.60, Obesity 23.28 ± 11.64, and MetS 23.00 ± 9.16. In contrast, low (Control 0.41 ± 0.70, Obesity 0.42 ± 0.81, MetS 0.78 ± 1.77) and intermediate heteroplasmic mtSNVs (Control 0.64 ± 1.41, Obesity 0.19 ± 0.40, MetS 0.96 ± 1.75) are relatively rare in CD14+ cells (Fig. 2). No statistically significant differences in mtSNV amounts were found between groups (low heteroplasmic mtSNV: P = 0.611; intermediate heteroplasmic mtSNV: P = 0.409; homoplasmic mtSNV: P = 0.955) (analysis was performed by Kruskal–Wallis test with Dunn’s multiple comparisons test).
mtSNV prevalence and frequency did not differ between groups; homoplasmies predominated
The heteroplasmic mtSNV signature did not display any distinctive patterns across the study groups; most heteroplasmic mtSNVs were unique events and occurred in fewer than 5% of individuals in each group (Fig. 3a). However, some individual heteroplasmies were more common. For example, among low heteroplasmic mtSNVs, the m.310T>C mutation was relatively frequent, with a prevalence of 17.64% in the Control and MetS groups and 9.52% in the Obesity group. Among intermediate heteroplasmic mtSNVs, the most common heteroplasmy in the Control and MetS groups was m.16189T>C (14.7% and 18.75%, respectively), while in the Obesity group, the most common heteroplasmy was m.15204T>C (9.52%). Alternative homoplasmies were the most common mutations. The most common homoplasmies, present in more than 90% of study participants, were found across all study groups: m.263A>G, m.750A>G, m.8860A>G, m.15326A>G, m.1438A>G, and m.4769A>G.
The frequency of various mtSNVs was analyzed across the entire mtDNA and its loci (Fig. 3b). The D-loop locus was the most variable mtDNA region in all three groups (Table 3). The tRNA locus had the lowest heteroplasmy frequency in all three groups compared to other loci (Table 3). The frequency of mtSNVs, both across the mtDNA and within specific loci, did not differ significantly between the study groups. Thus, the mtDNA of CD14+ cells in both patients and healthy controls showed similar frequencies of mtSNV variants across the mtDNA and at different loci.
Table 3Statistical description and group comparison of mtSNV frequencies (%)
| Locus | Category of mtSNV | Group | Median (Q1; Q3) | P-value |
|---|
| Whole mtDNA | Low heteroplasmic | Control | 0 (0; 0.01) | 0.2392 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0.01) | |
| Intermediate heteroplasmic | Control | 0 (0; 0.01) | 0.0568 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0.01) | |
| Homoplasmic | Control | 0.17 (0.08; 0.2) | 0.9548 |
| | Obesity | 0.18 (0.08; 0.21) | |
| | MetS | 0.17 (0.08; 0.18) | |
| D-loop | Low heteroplasmic | Control | 0 (0; 0.07) | 0.1857 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0.09) | |
| Intermediate heteroplasmic | Control | 0 (0; 0) | 0.0602 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0.09) | |
| Homoplasmic | Control | 0.49 (0.36; 0.69) | 0.9849 |
| | Obesity | 0.53 (0.36; 0.71) | |
| | MetS | 0.49 (0.36; 0.80) | |
| Protein coding | Low heteroplasmic | Control | 0 (0; 0) | 0.9682 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0) | |
| Intermediate heteroplasmic | Control | 0 (0; 0) | 0.2906 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0) | |
| Homoplasmic | Control | 0.12 (0.05; 0.16) | 0.9034 |
| | Obesity | 0.14 (0.04; 0.16) | |
| | MetS | 0.12 (0.05; 0.14) | |
| rRNA | Low heteroplasmic | Control | 0 (0; 0) | 0.1804 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0) | |
| Intermediate heteroplasmic | Control | 0 (0; 0) | 0.7231 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0) | |
| Homoplasmic | Control | 0.19 (0.15; 0.23) | 0.6223 |
| | Obesity | 0.19 (0.15; 0.19) | |
| | MetS | 0.19 (0.15; 0.19) | |
| tRNA | Low heteroplasmic | Control | 0 (0; 0) | 0.4234 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0) | |
| Intermediate heteroplasmic | Control | 0 (0; 0) | 0.4234 |
| | Obesity | 0 (0; 0) | |
| | MetS | 0 (0; 0) | |
| Homoplasmic | Control | 0.07 (0; 0.13) | 0.2139 |
| | Obesity | 0 (0; 0.07) | |
| | MetS | 0.07 (0; 0.07) | |
Associations between mtDNA variability and cardiometabolic parameters
Weak correlations were found between overall mtDNA frequencies, as well as at individual loci, and cardiometabolic parameters (Fig. 4a and b) after correction for confounders (sex, age, BMI) using the partial Spearman correlation test. The frequency of low heteroplasmies was associated with fasting blood glucose (FBG), blood urea nitrogen, low-density lipoproteins (LDL), and aspartate aminotransferase. For example, low heteroplasmy of the protein coding locus was correlated with LDL (r = −0.258; 95% CI −0.45 – −0.043), and D-loop low heteroplasmy correlated with FBG (r = 0.275; 95% CI 0.062–0.464).
Intermediate heteroplasmic mtSNV frequencies in the protein coding locus were also correlated with blood urea nitrogen (r = 0.297, 95% CI 0.08605–0.4828) and the degree of vascular stenosis (r = 0.396; 95% CI 0.067–0.647). Intermediate heteroplasmic mtSNV frequencies in the rRNA locus were negatively correlated with insulin levels (r = −0.228; 95% CI −0.424 – −0.019).
The highest number of correlations was found between homoplasmic mtSNVs and biochemical parameters. Homoplasmic mtSNV frequencies correlated with levels of alkaline phosphatase, total and indirect bilirubin, blood urea nitrogen, creatinine, and patient age (although the correlation with age was observed only for homoplasmic mtSNVs in the D-loop region).
Impact of mtDNA variability on the phenotypic parameters of CD14+ cells
To assess the impact of mtDNA variability on cellular phenotypic parameters, a multivariate PERMANOVA was conducted (all results presented in Supplementary Table 4). After adjusting for multiple comparisons, only one significant association was identified in human monocytes: intermediate heteroplasmic mtSNVs across the entire mtDNA were associated with the cytokine secretion stimulation index (FC): R2 = 0.1564 and P = 0.0030. The PERMANOVA results indicate that intermediate heteroplasmic mtSNVs in monocyte mtDNA account for approximately 15% of the secretory signature of these cells (R2 = 0.1564, P = 0.0030).
Discussion
Inspired by the approach of Calabrese et al.,15 we conducted a similar analysis, focusing on a single cell type that plays a key role in the pathogenesis of MetS components: blood monocytes. We considered not only individual donor diagnoses but also the cardiometabolic environment in which these monocytes reside. Therefore, next-generation sequencing was performed on 87 mtDNA samples from blood monocytes of MetS patients, obese patients, and healthy donors.
Monocyte mtDNA in our study showed weak mtDNA variability; the most common mutations were alternative homoplasmies, while heteroplasmies at both low and intermediate levels were almost absent. The distribution patterns of mtSNVs at mtDNA loci did not differ between patients and healthy controls. However, we identified some notable heteroplasmic hotspots present in more than 5% of the group: low heteroplasmic mtSNV m.310T>C, intermediate heteroplasmic mtSNVs m.16189T>C and m.15204T>C. None of these mtSNVs had any definitive pathological association in our study. However, for example, the m.16189T>C position in the D-loop has been repeatedly shown to be associated with the development of IR.21 In addition to heteroplasmic positions, we identified common alternative homoplasmies, some of which are foundational mutations of mitochondrial haplogroups. For example, according to the Mitomap database, mutations m.263A>G, m.750A>G, m.8860A>G, m.15326A>G, m.1438A>G, and m.4769A>G occur in at least 80% of lineages L, M, and N.22
Thus, at the genetic level, we did not find significant differences in mtSNVs between patients and healthy donors; however, the question of the influence of mtDNA variability on phenotypic indicators remains open.
We studied the associations between mtDNA variability and cardiometabolic parameters in donors. The protein-coding locus showed the most correlations, which is not surprising because it determines the main functions of mitochondria. Interestingly, a recent study on the prevalence of heteroplasmic mutations in various human tissues found that heteroplasmic mutations in the protein-coding locus occur most frequently in blood cells compared to other tissues.23
The rRNA locus had the next highest number of correlations. Only intermediate heteroplasmic mtSNVs correlated with biochemical parameters. This may indicate a dose-dependent effect of heteroplasmies in the rRNA locus: the higher the heteroplasmy level, the more significant its impact on cell phenotype.
One of the most interesting findings was a positive association between stenosis and intermediate heteroplasmies in the protein-coding locus: patients with a higher frequency of heteroplasmies in monocytes tended to develop stenosis. In contrast, low heteroplasmic mtSNVs in the protein-coding locus were negatively correlated with LDL levels. We hypothesize that mitochondrial heteroplasmy may result from active LDL uptake by blood monocytes and the development of vascular stenosis, rather than cause it.
Based on the obtained associations, we propose the following hypothesis: chronic LDL uptake by monocytes leads to mitochondrial dysfunction. Mitochondria cannot manage the large influx of substrate for oxidative phosphorylation, resulting in decreased mitochondrial membrane potential and the development of oxidative stress.24 Additionally, LDL accumulation in monocytes leads to the formation of foamy macrophages, a key step in atherosclerotic plaque formation.25 Thus, monocytes may help reduce blood LDL levels in patients with MetS, but LDL accumulation within monocytes can directly increase mtDNA variability, particularly the accumulation of heteroplasmy.
The D-loop was the only locus where the frequency of low heteroplasmic mtSNVs was associated with patients’ FBG levels. Negative correlations were also found between the frequency of intermediate heteroplasmic mtSNVs and insulin at the rRNA locus. This is an intriguing finding, as previous studies have shown that mutations in the tRNA locus (so-called tRNA modopathy), rather than rRNA loci, play a significant role in the development of IR, primarily by disrupting insulin secretion in beta cells.21,26 Among our 87 donors, none had tRNA mutations classically associated with IR.21
The role of rRNA mutations in insulin sensing has not been studied as extensively. Mouse strains with alternative mtDNAs are more prone to developing the MetS phenotype, which researchers attribute to alternative homoplasmies of mitochondrial rRNAs.27 Nevertheless, it seems confusing—how can monocyte mtDNA be associated with insulin sensitivity? In fact, monocytes have insulin receptors,28 and insulin signaling appears to actively influence their secretory potential. For example, short-term insulin stimulation of monocytes from healthy donors led to increased reactive oxygen species and IL-6 production, while chronic insulin exposure increased tumor necrosis factor production.29 Researchers also reported that IR affects the metabolic status of monocytes: monocytes from individuals with IR were relatively more reliant on oxidative metabolism.30 The observed associations between heteroplasmic mutations in monocyte mtDNA and blood insulin levels may be due to impaired insulin signaling in these cells, but this hypothesis remains to be confirmed.
While heteroplasmies correlated with “obvious” biochemical parameters in MetS, alternative homoplasmies were associated with less obvious parameters—alkaline phosphatase, creatinine, and bilirubin. As mentioned above, many of the alternative homoplasmies we identified are canonical mutations for specific haplogroups. Because haplogroups are an invariable parameter across different cell types, we may not be observing a monocyte-specific association between homoplasmies and biochemical parameters, but rather a likely organism-wide pattern. With further study involving mtDNA haplogroups (in this iteration of the study, our sample size is insufficient), we should be able to more fully describe the pathophysiological mechanisms driving these correlations.
Using a multivariate PERMANOVA analysis of variance, we found an association between intermediate heteroplasmic mtSNVs across the entire mtDNA and the cytokine stimulation index. This suggests that mtDNA heteroplasmy may influence monocyte secretory patterns, or conversely, that monocyte secretory characteristics may promote the accumulation of heteroplasmic mutations. However, it is worth noting that we did not find an association between the expression of genes of various groups (mitochondrial biogenesis, mitochondrial uncoupling, the oxidative stress system, and NF-kB signaling genes) and mtDNA variability. MtDNA variability may be associated with the expression of other signaling genes, such as cGAS-STING.31 Moreover, mtDNA variability may influence the formation of cytosolic mitochondrial RNA and DNA, which was not examined in this experiment.
To our knowledge, the role of mtDNA variability in cytokine secretion and immunobiogenesis in monocytes has not been widely studied. Most research has focused on mtDNA haplogroups rather than heteroplasmic mutations. In vitro experiments have shown that certain haplogroups are associated with increased gene expression of specific cytokines32,33; however, these findings have not been confirmed in primary cell cultures. A recent article examined the impact of an intermediate heteroplasmic mutation (m.5019A>G) with mtSNV levels between 70% and 87% in macrophages on cytokine secretion and phenotypic profile.34 The researchers showed that this heteroplasmic mutation causes mitochondrial dysfunction and alters the secretory profile of macrophages, resulting in increased type I interferon signaling and selective suppression of IL-1β and COX-2, as well as a slight decrease in IL-6 release. Notably, it is the intermediate heteroplasmy that leads to these phenotypic changes, which may be consistent with the results of our study. Certainly, these hypotheses require experimental confirmation and raise new questions: how exactly do mtDNA mutations regulate cytokine secretion? Through what processes and signaling pathways does this occur, and does it occur only in pathological or also in normal conditions?
The main results of the study are summarized in Figure 5.
The main limitation of the study is sample size, which is smaller than what is typical in population-based studies examining associations between genetic signatures and diseases. However, the goal was to holistically examine the place and role of genetic signatures across a broad spectrum of cardiometabolic parameters. Therefore, the small sample size is offset by the study’s comprehensive methodology, which makes it unique and relevant. It is important to note that our study demonstrates associations rather than causation and does not account for confounding factors such as diet and exercise. Another limitation is that, when analyzing mtDNA homoplasmic mutations, we did not clearly distinguish between haplogroup-related mutations and those that are functional or pathogenic. Most of the homoplasmic mutations we identified are haplogroup-related.