Variant load of mitochondrial DNA in single human mesenchymal stem cells

Abstract

Heteroplasmic mitochondrial DNA (mtDNA) variants accumulate as humans age, particularly in the stem-cell compartments, and are an important contributor to age-related disease. Mitochondrial dysfunction has been observed in osteoporosis and somatic mtDNA pathogenic variants have been observed in animal models of osteoporosis. However, this has never been assessed in the relevant human tissue. Mesenchymal stem cells (MSCs) are the progenitors to many cells of the musculoskeletal system and are critical to skeletal tissues and bone vitality. Investigating mtDNA in MSCs could provide novel insights into the role of mitochondrial dysfunction in osteoporosis. To determine if this is possible, we investigated the landscape of somatic mtDNA variation in MSCs through a combination of fluorescence-activated cell sorting and single-cell next-generation sequencing. Our data show that somatic heteroplasmic variants are present in individual patient-derived MSCs, can reach high heteroplasmic fractions and have the potential to be pathogenic. The identification of somatic heteroplasmic variants in MSCs of patients highlights the potential for mitochondrial dysfunction to contribute to the pathogenesis of osteoporosis.

Keywords: Bone disease; Mesenchymal stem cells; Osteoporosis; Somatic mtDNA variation.

Fig. 1

Histogram of somatic heteroplasmic variation in 13 patient samples. Frequency histogram of heteroplasmy fraction (HF) of all somatic heteroplasmic variants (HF > 0.00/< 1.00, n = 7849 variants) identified in 99 cells from 13 samples. The red line indicates a HF of 0.02. The vast majority of somatic heteroplasmic variation (n = 5,372 or 68.4%) had a HF < 0.02. 2,477 variants (31.6%) were taken forward for further analysis (SDataset 1).

Fig. 2

Somatic heteroplasmic count and heteroplasmic fraction distributions are similar across all 13 patient samples. Boxplots and strip charts showing the count of somatic heteroplasmies per cell per sample (upper) and the heteroplasmic fraction (lower, HF > 0.02/< 0.98) of each somatic heteroplasmic variant detected in each single-MSC grouped by sample. We observed no significant difference in heteroplasmy counts or heteroplasmy fraction distribution between samples (for each, ANOVA p > 0.05 and Dunnett’s p > 0.05 with patient 1 as reference). Boxes denoted 25th and 75th percentile, whiskers indicate 95^th percentile and the horizontal line indicates the median.

Fig. 3

Somatic heteroplasmic count and heteroplasmic fractions stratified by locus. Strip charts showing count of heteroplasmies (HF > 0.02/< 0.98) per cell per sample after normalisation by locus length (total region bp lengths: D-loop = 1124, rRNA = 2511, tRNA = 1486 and Protein Coding = 11,382) (left) and the heteroplasmic fraction (HF > 0.02/< 0.98) of somatic heteroplasmies (right) grouped by locus type. As expected, there is a higher frequency of heteroplasmic variants in the D-Loop compared to rRNA, tRNA and protein-coding regions (ANOVA p = 8.2 × 10^–5, Dunnett’s p, D-Loop v rRNA = 0.002, D-Loop v tRNA = 0.011, and D-Loop v protein coding = 0.001). There was no significant difference in overall HF between loci (ANOVA p = 0.115), However, rRNA HFs appear higher when compared to D-Loop as a reference (Dunnett’s p, D-Loop v rRNA = 0.038). tRNA and protein HFs distributions were not significantly different (Dunnett’s p, D-Loop v tRNA = 0.203 and D-Loop v protein coding = 0.082).

Fig. 4

Non Synonomous and synonomous somatic heteroplasmic count and heteroplasmic fractions stratified by locus. Upper Left: Strip chart showing the heteroplasmy (HF > 0.02/< 0.98) counts per cell per sample after normalisation by locus length (total region bp lengths: D-loop = 1124, rRNA = 2511, tRNA = 1486 and Protein Coding = 11,382). There was a significant increase in non-synonymous variant counts across all 13 individuals (Mann Whitney U, p < 0.001). Upper Right: Strip chart showing the distribution of protein-coding variant heteroplasmy fractions (HF > 0.02/< 0.98) after stratification into non-synonymous and synonymous variants. We observed no significant difference in HF between non-synonymous and synonymous variants (Mann Whitney U, p > 0.05). Lower left: Histogram of MutPred and APOGEE 2 scores showing a significant skew towards pathogenicity (MutPred score > 0.50, Kurtosis = 3.17, Shapiro-Wilks Normality Test p < 0.001 and APOGEE 2 score > 0.38, Kurtosis = 1.57, Shapiro-Wilks Normality Test p < 0.001). Dotted red lines indicate MutPred and APOGEE thresholds (0.5 and 0.38 respectively). Lower right: Boxplots of MSC somatic heteroplasmic MutPred and APOGEE 2 scores indicating that the majority of samples harbour cells with potentially pathogenic variants (MutPred pathogenic = 949 variants or 82% and APOGEE 2 = pathogenic 593 or 51%).