Characterization of the Mitochondrial Genome Landscape in MSCs, iPSCs and iMSCs from Osteoarthritis Patients and Healthy Donors

1. Introduction

Osteoarthritis (OA), the most prevalent joint disease worldwide, is predominately marked by synovial inflammation, subchondral bone remodeling and cartilage matrix degradation. Among OA risk factors, age over 50 years has been closely linked to OA pathogenesis contributing to inflammaging, cellular senescence, oxidative stress and impaired mitochondrial function [1]. As mitochondrial dysfunction is a hallmark of aging, the role of mitochondria in OA has attracted increased attention [2,3].

Mitochondria are multifunctional organelles that maintain cellular homeostasis by integrating bioenergetic processes with redox regulation and broader metabolic networks [4]. In addition to their dependence on nuclear-encoded proteins, mitochondria harbor their own (~16.5 kb) circular genome, rendering mitochondrial function on tight coordination between nuclear and mitochondrial genetic systems. Importantly, the functional consequences of heteroplasmic mutations are typically threshold-dependent with biochemical and cellular defects becoming evident only when mutant mtDNA exceeds a critical level [5].

In OA, mitochondrial dysfunction has been linked to alterations in the mitochondrial genome itself. Increased mtDNA damage, accumulation of somatic mtDNA variants and elevated heteroplasmy levels have been reported in OA-associated cell populations, such as chondrocytes and synoviocytes [3,6] suggesting that mtDNA instability may actively contribute to disease development and progression. Furthermore, oxidative stress-induced mtDNA alterations together with limited mtDNA repair capacity and clonal expansion of deleterious variants [7] are thought to exacerbate mitochondrial dysfunction, thereby establishing a pattern that promotes cartilage degeneration and impaired joint homeostasis.

Besides chondrocytes, other cells within the OA joint, as aged mesenchymal stromal cells (MSCs) display altered energy production, reduced oxidative phosphorylation efficiency and increased production of reactive oxygen species (ROS) contributing to cellular stress, phenotypic instability and reduced regenerative capacity [8]. Increasing evidence points that regulation of mitochondria dynamics and function is imperative for MSCs differentiation. Cell fate transitions, such as somatic cell reprogramming are accompanied by extensive metabolic and mitochondrial remodeling. During conversion of MSCs into induced pluripotent stem cells (iPSCs), a pronounced metabolic shift from oxidative phosphorylation toward glycolysis has been consistently reported, reflecting the bioenergetic requirements of pluripotency [9,10]. Re-differentiation of iPSCs into induced mesenchymal stromal cells (iMSCs) partially restores the mitochondrial phenotypic characteristics of primary MSCs. Although iMSCs regain aspects of mesenchymal stem cells’ identity and increase their reliance on oxidative metabolism, several studies have reported persistent differences in mitochondrial structure [11], bioenergetic function [12] and metabolic flexibility when compared to tissue-derived MSCs. Moreover, it has been reported that reprogramming to pluripotency is associated with substantial changes in mtDNA variant composition and heteroplasmy levels, including selective loss of pre-existing somatic mutations and emergence or expansion of new mtDNA variants [9]. The above observations support the existence of mtDNA selection mechanisms operating during cellular reprogramming [13] which may act to favor mitochondrial genomes compatible with the pluripotent state.

Until now, most studies in OA cell populations have primarily examined the global mitochondrial mutation burden without resolving how mtDNA variants are distributed across individual mitochondrial genes or respiratory chain complexes [14], nor whether cellular reprogramming effectively resets OA-associated mitochondrial genomic alterations [9,15].

In the present study, we systematically characterized mtDNA heteroplasmy and variant distribution across bone marrow derived-MSCs, iPSCs and iMSCs obtained from OA patients and young healthy donors. By integrating deep mitochondrial genome sequencing with transcriptomic and miRNome data, we aimed to investigate mitochondrial remodeling during cellular reprogramming and assess mitochondrial genome interaction with nuclear-encoded mitochondrial genes in relation to OA environment.

2. Materials and Methods

2.1. Bone Marrow Samples

Bone marrow samples were collected from 3 osteoarthritis (OA) patients and 2 healthy donors. The OA cohort was comprised of two females and one male with ages 69, 76 and 79 years, respectively. Bone marrow aspirates were obtained from the proximal metaphysis of the femur immediately following femoral neck osteotomy during total hip arthroplasty at the Orthopaedics Department of the University Hospital of Larissa, Greece. Healthy donors, one female and one male, ages 6 and 13 years, were pediatric cell-transplant volunteers at the Bone Marrow Transplantation Unit of “Aghia Sophia” Children’s Hospital in Athens, Greece and bone marrow aspirates were collected from the medial iliac crest. The study protocol was approved by the Ethics Committees of both “Aghia Sophia” Children’s Hospital and University Hospital of Larissa and complied with the principles of the Declaration of Helsinki (1975). Informed consent was obtained from all OA patients and from the parents or legal guardians of the underaged healthy donors.

2.2. BM-MSCs Isolation, iPSC Generation and Differentiation into iMSCs

Isolation of BM-MSCs, reprogramming into iPSCs, and subsequent differentiation into induced MSCs (iMSCs) followed established protocols previously described by us [16,17]. Briefly, BM-MSCs from OA patients (MSC_OA) and healthy donors (MSC_N) were expanded in DMEM supplemented with 10% FBS, 1% Pen/Strep, HEPES, L-glutamine, ascorbic acid, sodium pyruvate and L-glucose under standard incubation conditions (37 °C, 5% CO₂). Adherent mesenchymal-like cells appeared within the first days of culture and were sub-cultured upon confluence using trypsin at a 1:3 split ratio. Reprogramming of MSCs into iPSCs was performed using a modified synthetic mRNA-based approach, as previously reported [16]. In brief, MSCs were seeded onto irradiated NuFF feeder cells and repeatedly transfected with the five canonical reprogramming factors (OCT4, SOX2, KLF4, c-MYC and LIN28) in the presence of B18R. Emerging colonies with pluripotent morphology were manually selected and expanded in mTeSR1 medium on Matrigel-coated surfaces. For mesenchymal induction, iPSCs were subjected to embryoid body formation followed by adherence to gelatin-coated flasks in MSC differentiation medium, as previously described [16,17]. Resulting induced mesenchymal stem cells from OA patients (iMSC_OA) and healthy donors (iMSC_N) lines were subsequently expanded under the same culture conditions used for primary MSCs. Characterization of iMSCs included flow cytometric profiling of canonical MSC markers (CD90, CD44, CD105, CD73, CD29) and exclusion of hematopoietic markers (CD34, HLA-DR). Cell counts and viability were assessed using 7-aminoactinomycin D (7-AAD). Trilineage differentiation capacity (osteogenic, adipogenic and chondrogenic) was confirmed according to the manufacturers’ protocols using lineage-specific differentiation media followed by histochemical staining [16].

2.3. DNA Extraction and Mitochondrial DNA Amplification

Genomic DNA was extracted from MSC_OA (n = 3), MSC_N (n = 2), iPSC_OA (n = 3), iPSC_N (n= 2), iMSC_OA (n = 3) and iMSC_N (n = 2) samples using a Qiagen DNA extraction kit (Qiagen, Hilden, Germany) according to manufacturer’s instructions. DNA concentration and purity were assessed using NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and Qubit 4.0 fluorometer (Thermo Fisher, Waltham, MA, USA), while genomic integrity was evaluated by electrophoresis on 1% agarose gels. To enrich and linearize the mitochondrial genome, long-range PCR was performed to amplify two overlapping fragments spanning the entire mitochondrial sequence (rCRS, GenBank ID NC_012920). Amplicons were generated using primer sets previously described [18]: Set 1 (5042f–1424r, 12.96 kb) and Set 2 (528f–5789r, 8.7 kb). Primer specificity was validated by the absence of amplification in Rho0 control DNA, confirming that nuclear mitochondrial pseudogenes (NUMTs) were not co-amplified. PCR reactions were performed using LongAmp Taq DNA Polymerase (New England Biolabs, Ipswich, Massachusetts, USA) in 25 μL reaction volumes containing 50 ng of genomic DNA. Amplified products were assessed on 0.8% agarose gels alongside positive and negative controls and subsequently quantified using a Qubit 4.0 fluorometer. Purification of PCR products was carried out using AMPure XP beads (Beckman Coulter, Brea, CA, USA). Following purification, amplicons were pooled in defined equimolar ratios (0.35 × Set 2 and 0.65 × Set 1) to ensure proportional representation of each fragment. The pooled mixture was re-quantified and subjected to a second purification step using the QIAquick PCR Purification Kit (Qiagen, Hilden, Germany) prior to downstream library preparation and sequencing.

2.4. Mitochondrial DNA Deep-Sequencing

Pooled amplicons were tagmented, amplified, purified, normalized and subsequently combined into a single Illumina library using the Nextera XT DNA Sample Preparation Kit (Illumina, CA, USA), according to manufacturer’s instructions. For library preparation, 1 ng of DNA per sample was used. Resulting libraries were quantified using Qubit 4.0 Fluorometer (Thermo Fisher, Waltham, MA, USA) and evaluated for fragment size distribution and overall quality using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Libraries were then normalized to equimolar concentrations and pooled into a multiplexed library suitable for high-throughput sequencing. Sequencing was performed on an Illumina MiSeq platform employing the MiSeq Reagent Kit v3.0 (Illumina, CA, USA) to generate paired-end reads of 2 × 150 bp. Sequencing runs were performed to a minimum depth of approximately 6000Χ per sample, enabling robust detection of low-frequency mitochondrial variants and accurate estimation of heteroplasmy levels.

2.5. Variant Calling and Bioinformatic Analysis

Raw sequencing reads were subjected to quality control using FastQC (v0.11.8) and subsequently aligned to the human mitochondrial reference genome (NC_012920.1) using the bioinformatic tool BWA (v0.7.17). The resulting alignment files were sorted and indexed with Samtools (v1.12), and PCR duplicates were identified and marked with Picard (v2.2.4). Mitochondrial single nucleotide variants (mtSNVs) were detected using Mitoverse online platform (mtDNA-Server v2 Beta (rc12), v2.0.0), which incorporates Haplogrep and Haplocheck packages for variant calling, haplogroup assignment and evaluation of potential sample contamination. Heteroplasmy fraction (HF) was calculated, as the ratio of variant allele depth to reference allele depth. Variants with HF>0.98 were classified as homoplasmic, whereas HF values >0.02 and <0.98 were considered heteroplasmic. Heteroplasmic burden was defined as the total number of heteroplasmic variants detected per sample. To ensure high-confidence calls we applied the following filtering criteria: (1) variant allele fraction (AF) >2%; (2) mean mtDNA depth of coverage >1500X to enable detection of low-level heteroplasmy; (3) retention of single-nucleotide variants only; (4) exclusion of variants exhibiting strand bias and (5) removal of variants located in known low-complexity mtDNA regions (66–71 bp, 300–316 bp, 513–525 bp, 3106–3107 bp, 12,418–12,425 bp and 16,181–16,194 bp). To exclude sequencing-related artifacts, mean mitochondrial read depth was compared between OA and healthy donors. No statistically significant difference was detected (unpaired t-test) indicating that the observed heteroplasmy differences were unlikely to be driven by sequencing-depth variation. Putative pathogenicity of mtSNVs was evaluated using the MutPred algorithm, and variants scoring >0.70 were classified as potentially pathogenic in accordance with the MITOMAP Human Mitochondrial Genome Database.

2.6. Transcriptomic Data Analysis

Transcriptomic (mRNA and lncRNA) and miRNome datasets, previously generated in our laboratory [17] were reanalyzed for the purposes of this study. Differential expression analysis was performed on filtered datasets using DESeq2-normalized values [19]. The false discovery rate (FDR) was controlled using the Benjamini–Hochberg correction. Genes were considered differentially expressed when the adjusted p-value was ≤0.05 and the fold-change (FC) was >1.5 or <-1.5. Venn diagrams and heatmaps were generated using the R/Bioconductor packages VennDiagram and ComplexHeatmap, respectively. Nuclear-encoded mitochondrial genes were identified using the Human MitoCarta3.0 database [20]. Regulatory interactions between differentially expressed miRNAs and their predicted target genes and lncRNAs were visualized using the network package in R and the Cytoscape platform (v3.10.1).

2.7. Statistical Analysis

Statistical analyses were performed to evaluate the effects of donor status (OA vs healthy) and cell type (MSCs, iPSCs, iMSCs) on mitochondrial variant burden and related parameters. For within-donor comparisons across the three cell states, non-parametric repeated-measures tests were used. Specifically, one-way ANOVA and/or Friedman’s test was applied to compare total heteroplasmic variant counts and hotspot gene variant burden across MSCs, iPSCs and iMSCs derived from the same donors, followed by Wilcoxon signed-rank post-hoc tests with Benjamini–Hochberg correction for multiple pairwise comparisons, where appropriate. For two-group comparisons between OA and healthy donors within the same cell type, Mann-Whitney U test was used. Mean mitochondrial sequencing depth was compared between OA and healthy samples using either unpaired t-test or Mann–Whitney U test, depending on the distribution of the data, to exclude depth-related bias in heteroplasmy estimates. Unless otherwise stated, statistical significance was defined as p<0.05. All statistical analyses were conducted in R (v3.10) using standard statistical packages.

3. Results

3.1. Reprogramming to iMSCs Reduces the Heteroplasmic Variant Burden of OA-Derived MSCs

Quantitative assessment of mitochondrial single nucleotide variants (SNVs) across the three cell populations (MSCs, iPSCs, iMSCs) derived from osteoarthritis (OA) patients and healthy donors revealed clear differences in heteroplasmic variant burden among them. Specifically, we found that primary bone marrow MSC_OA exhibited markedly higher number of heteroplasmic variants (Heteroplasmic Fraction, HF>2%) compared to MSC_N. On average, MSC_OA harbored 198-335 heteroplasmic SNVs per sample, whereas MSC_N displayed only 7-15 variants (Mann-Whitney U test, p=0.015, Figure 1A and Figure 1B). Following cellular reprogramming and subsequent differentiation, progressive reduction in heteroplasmic load was detected in the OA group. Specifically, iPSC_OA showed substantial decrease in the number of variants, while the corresponding iMSC_OA exhibited a slight further decline (one-way ANOVA, p= 0.012, Figure 1A). Notably, within the OA cohort, inter-donor variability was evident primarily among MSCs (one-way ANOVA, p= 0.042, Figure 1A), indicating that mitochondrial heterogeneity in OA may be donor dependent. In contrast, in cell populations derived from healthy donors, the heteroplasmic load remained consistently low and stable across all cellular states (MSCs iPSCs, iMSCs) with no significant differences among them (one-way ANOVA, p = 0.329) (Figure 1B).

In addition to the quantitative assessment of heteroplasmic variant burden and allele frequencies, we also examined the qualitative distribution of SNVs across the mitochondrial genome. We found that in primary MSC_OA, the majority (60.67%) of heteroplasmic variants were located within protein-coding regions, suggesting a potential direct effect on mitochondrial protein synthesis and OXPHOS function. The predominance of coding-region variants persisted, albeit at lower proportions, following reprogramming into iPSCs (45.78%) and subsequent differentiation into iMSCs (56%). Notably, we observed that at the iMSC stage no variants were detected within tRNA genes, while a relative increase in variants mapped to regulatory non-coding regions, particularly the control region (CR). The proportion of variants localized in rRNA genes remained relatively constant across all OA-derived cell types, representing approximately 25% of the total heteroplasmic variants (Figure 1C). A comparable localization pattern was observed in all cell populations derived from healthy donors, although the overall number of heteroplasmic variants was markedly lower. Specifically, we found that in MSC_N, approximately 60% of variants were found in coding regions, while very few mapped to regulatory sequences. During reprogramming and redifferentiation, as the total heteroplasmic load further declined, no tRNA-associated variants were detected (Figure 1D).

3.2. Reprogramming Eliminates Variant Hotspots in OXPHOS and Mitochondrial Ribosomal Genes in OA-Derived iMSCs

We next seek to investigate the heteroplasmic SNV accumulation per mitochondrial gene in the OA derived cell populations (MSCs, iPSCs, iMSCs). We observed that in MSC_OA certain mitochondrial genes emerged as distinct “hotspots” of SNV accumulation. Notably, MT-ND1, MT-ND2, MT-ND3 and MT-ND5 genes, which encode subunits of Complex I of the respiratory chain, MT-ATP8 encoding Complex V, an ATP synthase subunit and MT-CO1 encoding Complex IV, a cytochrome c oxidase subunit 1, harbored on average approximately 30 heteroplasmic variants each (Figure 2). Following reprogramming of MSCs into iPSCs and their subsequent differentiation into iMSCs, a marked reduction in variant load to less than five variants per gene was observed across the above mentioned 6 genes (Friedman’s test, MT-ND1, adjusted-p=0.027, MT-ND2, adjusted-p=0.032, MT_ND3, adjusted-p=0.0193, MT-ND5, adjusted-p=0.0441, MT-ATP8, adjusted-p=0.047, MT-CO1, adjusted-p=0.036). A similar pattern was evident in the ribosomal RNA genes MT-RNR1 and MT-RNR2, which encode the 12S and 16S rRNAs, respectively. In primary MSC_OA, MT-RNR1 and MT-RNR2 exhibited particularly high variant counts (approximately 40 SNVs), while following reprogramming and redifferentiation, nearly all rRNA variants were eliminated, approaching baseline levels comparable to those observed in healthy controls (Friedman’s test, MT-RNR1, adjusted-p=0.041, MT-RNR2, adjusted-p=0.046, Figure 2). In samples derived from healthy donors, the distribution of mitochondrial variants across the different cell types did not demonstrate any distinct pattern. Overall variant burden remained low and no recurrently altered regions or mutational “hotspots” were observed (Supplementary Figure 1)

3.3. Differentiation Selectively Eliminates Highly Pathogenic Mitochondrial Variants in OA-Derived iMSCs

To assess the potential functional impact of the heteroplasmic variants identified in MSC_OA, the protein-coding SNVs were stratified into synonymous and non-synonymous (Figure 3A). Across all cell types, differences in variant composition were observed between synonymous and non-synonymous substitutions. In both iPSCs and iMSCs, synonymous variants outnumbered non-synonymous variants, indicating a relative predominance of functionally neutral substitutions. Analysis of heteroplasmic fractions (HFs) revealed no significant differences between variant classes in iPSCs and iMSCs, although synonymous variants tended to exhibit higher mean HFs (Mann–Whitney U test, iPSCs p = 0.137; iMSCs p = 0.42). In contrast, within MSC_OA, non-synonymous variants were more abundant than synonymous variants (211 vs 145), reflecting a shift toward protein-altering mitochondrial variation in the diseased state (HFs between synonymous and non-synonymous variants, Mann–Whitney U test, p = 0.0431)

To further investigate the above observation, we analyzed the origin and persistence of each variant by comparing the complete SNV profile among all OA cell types (Figure 3B). Variants unique to each cellular population were subsequently assessed for the relationship between their heteroplasmic fraction and predicted pathogenic potential using the MutPred algorithm. We observed that among unique variants detected in MSC_OA (n=281, Supplementary Table 1), which were absent after reprogramming and redifferentiation, mutations with high MutPred scores (>0.6) indicative of elevated pathogenic potential, tended to exhibit relatively high heteroplasmic fractions. Importantly, in MSC_OA samples, pathogenicity was positively correlated (slope = 0.153) with heteroplasmic fraction implying that more deleterious variants tended to accumulate at higher frequencies (Figure 3C). Although most variants displayed low heteroplasmy (<5%), two notable exceptions, m.7913C and m.7821T were identified, which demonstrated relatively high heteroplasmy (~18%) and elevated MutPred score (~0.80).

Moreover, three additional variants, m.3668T, m.14749G, and m.4789T, displayed very high predicted pathogenicity (~0.95) despite low heteroplasmic fractions, suggesting potential functional impact even at subthreshold levels. In contrast, 15 unique variants (Supplementary Table 2) detected exclusively in iMSCs_OA were characterized by very low heteroplasmy (~2%) and minimal predicted pathogenicity, indicating limited biological significance. Exceptions were the variants m.9939A and m.1399C which showed marginal MutPred scores but remained at low heteroplasmy levels suggesting negligible functional effect (Supplementary Figure 2).

3.4. Regulation of Nuclear-Encoded Mitochondrial Genes Reveals Metabolic Reprogramming During the MSC_OA to iMSC_OA Transition

To identify nuclear-encoded genes functionally linked to mitochondrial activity, the previously identified by us differentially expressed genes (DEGs) in iMSCs [17] were cross-referenced with MitoCarta 3.0 database. The comparison revealed five genes, COX7A1, GPAT2, MAOB, SERAC1, and ECHDC2 that were directly associated with mitochondrial function (Figure 4A). These five genes exhibited markedly reduced expression in iMSCs_OA resembling that of iPSCs and MSCs derived from healthy donors (MSC_N) (Figure 4C).

We then constructed a lncRNA-miRNA-mRNA competing endogenous RNA (ceRNA) interaction network incorporating exclusively the identified nuclear genes and their interacting regulatory RNAs in OA-derived cell types (Figure 4B). Network analysis revealed complex interactions between 28 hub miRNAs, as let7b/i, miR-497-5p and 9 hub lncRNAs, as MEG3/8, SNHG14 displaying inverse expression patterns. We observed that the iMSCs-lncRNA profile formed an intermediate state between iPSCs and MSC_OA indicating partial reset of disease-associated post-transcriptional signatures following reprogramming.

4. Discussion

Recent evidence has revealed that mitochondrial dysfunction in aged cells, manifested as mtDNA damage, altered mitochondrial dynamics and metabolic imbalances, contributes to several diseases, as cancer, neurodegenerative, cardiovascular, liver and metabolic diseases [21,22]. A significant role of mitochondrial dysfunction has also been attributed to the onset and progression of osteoarthritis [3].

In the present study, using ultra-deep sequencing we evaluated single-nucleotide variants (SNVs) in MSCs, iPSCs, and iMSCs populations derived from healthy donors and OA patients. MSCs from OA patients (MSC_OA) revealed elevated heteroplasmy reflecting mitochondrial genomic instability associated with aging and chronic oxidative stress [23], conditions that are present in OA [24]. The observed increase in heteroplasmic SNVs counts indicates cumulative mitochondrial damage driven by oxidative load and the inflammatory microenvironment of OA. Prior studies have demonstrated that aged MSCs harbor increased mtDNA SNV burden, particularly in genes associated with Complex I and Complex IV contributing to reduced oxidative phosphorylation efficiency and impaired cellular energy output [11,25]. Additionally, we observed progressive reduction of heteroplasmic variants during reprogramming into iPSC_OA and iMSC_OA representing a mitochondrial bottleneck phenomenon, whereby reprogramming favors the selection of compatible mtDNA subpopulations leading to elimination of low-frequency variants and restoration of mitochondrial function. This observation is in accordance with recent studies reporting that iPSCs and iMSCs retain a remodeled mitochondrial genome with lower mutational burden accompanied by energetic recovery through re-activation of OXPHOS and improved ATP-production capacity [12,26]. Accordingly, the observed mitochondrial “rebalancing” in iMSC_OA may represent a rejuvenation-associated mechanism through which mitochondrial quality and energy output returns to near-physiological levels.

Regarding qualitative characteristics in the mitochondrial genome, we found that in MSC_OA certain mitochondrial genes as, MTND1/2/3/5MT-ATP8, and MT-CO1 encoding subunits of complexes I, V, and IV of the oxidative phosphorylation (OXPHOS) system emerged as distinct “hotspots” of heteroplasmic SNVs accumulation. The above novel finding is of particular interest, as these SNVs hotspots affect critical components of the OXPHOS machinery, potentially compromising mitochondrial energy production in MSCs derived from OA patients and contributing to cellular dysfunction in OA. Similar patterns have been reported in MSCs and chondrocytes isolated from OA tissues, where complex I dysfunction has been linked to increased reactive oxygen species (ROS) generation, disruption of mitochondrial membrane dynamics, and accelerated cellular senescence [27,28,29]. Following MSC_OA differentiation into iPSC_OA and subsequently into iMSC_OA, we observed elimination of variants localized to “SNVs hotspots”. The observed reduction of SNVs in MT-ND, MT-CO, and MT-ATP8 genes, as well as in MT-RNR1 and MT-RNR2 encoding the 12S and 16S rRNAs, is indicative of a partial “genetic reset” of mitochondrial genomic integrity in iPSCs and iMSCs. This trend suggests that iPSCs reprogramming and iMSCs differentiation exert a purifying effect on the mitochondrial genome, leading to partial restoration of heteroplasmic stability. This supports the hypothesis that such mutations may exert negative selective pressure during the acquisition of pluripotency and lineage-specific differentiation [30], potentially hindering mitochondrial performance and cellular reprogramming efficiency.

Interestingly, despite the increased number of SNVs observed in primary MSC_OA, the heteroplasmic fraction of most variants remained low (2–5%), likely indicating a polyclonal mitochondrial population, in which multiple low-frequency mtDNA genotypes coexist without a single dominant genome being positively selected [31]. In all reprogrammed cell states, iPSCs and iMSCs of healthy donors, we observed that synonymous variants exhibited higher heteroplasmic fractions consistent with neutral drift, rather than selective pressure. In contrast, in MSC_OA, an enrichment of non-synonymous variants was observed, which allows potentially deleterious substitutions to persist and expand within the mitochondrial population. This pattern, which is in line with the chronic oxidative and inflammatory environment of OA, may reflect disrupted mitochondrial surveillance mechanisms, including mitophagy, representing a cellular response to chronic damage and aging.

Based on the above findings, we then proceeded to identify non-synonymous mtDNA variants unique to MSC_OA and evaluate their potential functional impact. We identified a subset of variants (n=281) exclusive to MSC_OA that were no longer detectable following reprogramming. Notably, variants with higher predicted pathogenicity, as indicated by elevated MutPred scores (>0.6) were associated with relatively increased heteroplasmic fractions. Furthermore, the positive correlation observed between heteroplasmic fraction and predicted pathogenic potential among MSC_OA-specific variants suggests that potentially deleterious substitutions can accumulate to appreciable levels in the OA cellular state, supporting a defect of purifying selection mechanisms prior to reprogramming [32,33]. Notably, the identification of two rare, pathogenic according to MITOMAP [34] and the literature [35,36], MT-CO2 variants, m.7821C>T and m.7913T>C with elevated MutPred scores and relatively high heteroplasmic fraction underscores the potential for direct impairment of complex IV activity with downstream consequences for ATP production and ROS generation. However, the observed persistence in MSC_OA cells of the additional highly pathogenic variants, m.3668T, m.14749G, and m.4789T at low heteroplasmic fractions suggests that even low-threshold mtDNA alterations may contribute cumulatively to mitochondrial dysfunction in OA-derived MSCs [11]. Importantly, the complete absence of high-heteroplasmy pathogenic variants in iPSC_OA and iMSC_OA is consistent with a model of purifying selection during reprogramming and redifferentiation processes.

Taking into consideration that mitochondrial function is not only regulated by the mitochondrial genome, but nuclear-encoded genes, also, play critical roles in mitochondrial functionality, we proceeded, for the first time to our knowledge, to examine the behavior of the nuclear-encoded mitochondrial-associated genes. Cross-reference of the previously found by us differentially expressed genes (DEGs) in iMSCs during the differentiation process from MSCs_OA to iMSCs [17] with nuclear genes encoding mitochondrial proteins from the MitoCarta database [20], revealed 5 genes, COX7A1, GPAT2, MAOB, SERAC1, and ECHDC2, that were associated with mitochondrial function. All 5 nuclear genes exhibited reduced expression in iMSCs adopting an expression profile that closely resembled that of iPSCs and MSCs derived from healthy donors rather than MSC_OA samples. The observed reduced expression of COX7A1 and MAOB in iMSCs may indicate an adaptive down-tuning of mitochondrial oxidative activity to avoid excess ROS generation, profile compatible with post-rejuvenation metabolic stabilization state [37,38]. COX7A1, a subunit of Complex IV, has been linked to oxidative stress responses in aged or inflammatory cells [39,40], while MAOB inhibition has been associated with mitochondrial protection and improved stem cell energetic output [41,42]. In parallel, genes such as SERAC1, GPAT2, and ECHDC2, which regulate the phospholipid recycling and the architecture of the inner mitochondrial membrane system [43], have established roles in mitochondrial membrane integrity and OXPHOS efficiency. SERAC1 is essential for remodeling phosphatidylglycerol species and cardiolipin, a key structural component of the mitochondrial membrane and its dysfunction is associated with impaired mitochondrial biogenesis [44]. GPAT2, a regulator of lipid homeostasis, has been reported to have reduced expression in primary MSC populations [45]. The observed suppression of GPAT2 in iMSCs may reflect metabolic rebalancing following removal of OA-associated oxidative stress. In that regard, the combined downregulation of the nuclear-encoded mitochondrial-related genes may, therefore, indicate an adjustment in lipid-biosynthetic mechanisms to restore mitochondrial membrane stability and reduce metabolic fragility.

Finally, the constructed lncRNA-miRNA-mRNA interaction networks for COX7A1, GPAT2, MAOB, SERAC1, and ECHDC2 revealed a complex post-transcriptional regulatory framework. Multiple lncRNAs, such as MEG3 and SNHG14 may function as competing endogenous RNA (ceRNA) “sponges” targeting miRNAs [46,47], including hsa-miR-145, hsa-miR-195 and hsa-let-7b family members, ultimately modulating their gene targets. Integration of these multi-omic layers suggests that post-transcriptional RNA interactions contribute to mitochondrial phenotypic adaptation in iMSCs promoting energetic efficiency and limiting metabolic stress [48]. The above findings indicate that the observed mitochondrial genome “purification” in iMSCs is not an isolated genomic event but occurs in parallel with a coordinated reprogramming of nuclear-encoded mitochondrial pathways. The deregulation of these genes reflects a multi-omic metabolic recalibration, whereby post-transcriptional regulatory networks actively shape mitochondrial homeostasis supporting the restoration of functional efficiency following reprogramming.

Collectively, our findings provide novel insights into the fidelity of mitochondrial genome remodeling during cellular reprogramming. Such remodeling likely underpins the metabolic and functional reactivation observed in iMSCs, potentially contributing to the restoration of mitochondrial energy homeostasis and regenerative capacity. The above have important implications for the use of reprogrammed cells in OA modeling and support their potential applicability in regenerative medicine approaches targeting mitochondrial dysfunction.