Mitochondrial DNA Haplogroup N9a Negatively Correlates with Incidence of Hepatocellular Carcinoma in Northern China.

Mitochondrial DNA (mtDNA) haplogroups are associated with various types of cancer; however, the molecular mechanisms by which mtDNA haplogroups affect primary hepatocellular carcinoma (HCC) are not known. In this study, we carried out a case-control study on 388 HCC patients and 511 geographically matched asymptomatic control subjects in northern China. We found that mtDNA haplogroup N9a and its diagnostic SNP, m.16257C > A, negatively correlated with the incidence of HCC in northern China (odds ratio [OR] 0.290, 95% confidence interval [CI] 0.123-0.685, p = 0.005), particularly in patients with infection of hepatitis B/C virus (HBV/HCV) (for haplogroup N9a: OR 0.213, 95% CI 0.077-0.590, p = 0.003; for m.16257C > A: OR 0.262, 95% CI 0.107-0.643, p = 0.003). However, mtDNA haplogroup N9a is not associated with clinical characteristics of HCC including serum alpha-fetoprotein (AFP) level and tumor size. In addition, cytoplasmic hybrid (cybrid) cells with N9a haplogroup (N9a10a and N9a1) had transcriptome profiles distinct from those with non-N9a (B5, D4, and D5) haplogroups. Gene set enrichment analysis (GSEA) showed that metabolic activity varied significantly between N9a and non-N9a haplogroups. Moreover, cells with haplogroup N9a negatively correlated with cell division and multiple liver cancer pathways compared with non-N9a cells. Although it is still unclear how N9a affects the aforementioned GSEA pathways, our data suggest that mtDNA haplogroup N9a is negatively correlated with the incidence and progression of HCC in northern China.

Introduction

In China, liver cancer is the second most common malignancy, with an estimated 360,000 new cases and 350,000 deaths per year. Hepatocellular carcinoma (HCC) is the most common subtype of primary liver cancer (>80%). Most HCC is caused by chronic liver injury due to hepatitis B/C virus (HBV/HCV) infection and alcohol abuse. Although the molecular mechanisms underlying the development and progression of HCC are poorly defined, functional studies suggest that both HBV/HCV infection and alcohol abuse are closely associated with mitochondrial dysregulation,3, 4 which has long been associated with the development of HCC.5, 6

Mitochondria generate the energy necessary to sustain all cellular activities. Human mitochondria contain approximately 1,500 proteins, of which 13 are encoded by mtDNA. mtDNA variations are associated with multiple human conditions such as aging, cancer, and rare mitochondrial diseases. Moreover, mtDNA haplogroups, specific mtDNA genetic variations defined by germline mtDNA mutations retained through evolution, are associated with cancer and type 2 diabetes, as well as neurodegenerative diseases, such as Leber hereditary optic neuropathy12, 13 and Alzheimer’s disease. Mechanistically, different mtDNA haplogroups can differentially affect mitochondrial performance. For example, haplogroups M and N have different mitochondrial matrix pH and intracellular calcium levels; haplogroup J is more transcriptionally efficient and displays a higher level of mtDNA replication than haplogroup H. Transcriptome analysis revealed significant differences in mRNA levels among cytoplasmic hybrid (cybrid) cells of different mtDNA haplogroups, and by using mouse models of various mtDNA haplogroups, A. Latorre-Pellicer et al. saw that liver function varied significantly among mice with different mtDNA genetic backgrounds.

Three recent studies, including participants from northern China (Hebei), eastern China (Shanghai), and southern China (Guangxi), demonstrated that distinct mtDNA haplogroups are associated with various HCC-associated factors.18, 19, 20 Notably, effects of mtDNA haplogroup defining SNPs on the prognosis of HCC were conflicted,19, 20 and the correlations between specific mtDNA haplogroup and the incidence of HCC were not ever investigated in these studies. Multiple factors may contribute to this unexpected finding. Different environmental conditions among participants may have differentially affected the role of mtDNA haplogroups in HCC.

Pitfalls including insufficient or inappropriate recruitment criteria of participants, improper choice of statistical analysis methods, and lack of replication cohort may have also cast doubt on the validity of the conclusion.21, 22 In addition, two of the studies included patients with small sample size in their analyses;18, 20 therefore, their findings may not reflect the potential association of mtDNA SNPs/haplogroups with HCC-associated factors. To date, only the study using eastern Chinese population adopted a considerable number of patients to analyze the relationships between mtDNA haplogroup and HCC. To verify and to further assess the possible contribution of mtDNA haplogroup on the prevalence of HCC, we conducted a cohort study with participants from northern China to reevaluate the potential correlation between mtDNA haplogroups and the incidence and severity of HCC.

Results

mtDNA Haplogroup N9a Is Associated with Decreased Risk of HCC

To test whether there is a population stratification between the patients and control subjects, which may tone down the confidence level of the study, we performed a principal-component analysis (PCA) for patients and control subjects based on mtDNA haplogroup frequencies as previously described. As shown in Figure 1, both patients and control subjects are close with each other, and were co-clustered with most of the northern Chinese population. This result suggested that there is no population stratification between the patients and control subjects.

Figure 1

The PCA of HCC Patients and Controls, as well as Other Reported Han Populations across China, Based on mtDNA Haplogroup Frequencies

HCC patients and controls are from Zhengzhou, Henan, are shown with yellow and green fonts, the previously reported populations from Henan are shown with red fonts, and the previously reported populations from northern and southern China are shown with blue and purple fonts.

While mtDNA haplogroup M was found to be associated with late onset of HCC in northern China with limited sample size, we first asked whether the frequency of haplogroup M associated with HCC in larger sample size. However, we found that haplogroup M is not associated with HCC with adjustment of age and sex (Table 1) and is not associated with HCC at late onset by using the Kaplan-Meier method (Figure 2). Furthermore, no significant difference was found in the prevalence of macro-haplogroups A, B, D, F, G, and M7–10 between patients and controls (Table 2). We did, however, find that the prevalence of haplogroup N9 was significantly lower in patients with HCC compared to controls, using multivariate logistic regression analysis with adjustment for age and sex (odds ratio [OR] 0.241, [95% confidence interval (CI) 0.110–0.531], p < 0.001) (Table 2). Further analysis of sub-haplogroups from ten macro-haplogroups revealed that sub-haplogroup N9a (OR 0.248, [95% CI 0.101–0.610], p = 0.002) but not sub-haplogroup Y of N9 contributes to the negative association between macro-haplogroup N9 and HCC (Table 2).

Table 1

Analysis of Association between Macro-Haplogroup M/N and HCC

Macro-Haplogroup	Patients (n = 388)	Controls (n = 511)	OR (95% CI)	p Value
M	229 (59.0)	278 (54.4)	1.154 (0.875–1.520)	0.310
N	159 (41.0)	233 (45.6)
Age			1.000 (0.989–1.011)	0.992
Gender			0.322 (0.238–0.436)	0.000∗

p < 0.025 (0.05/2), adjusted p value with Bonferroni correction while two haplogroups were studied. Values in parentheses are the percentage of samples. CI, confidence interval; OR, odds ratio.

Figure 2

The Association between Macro-Haplogroup M/N and Age at Onset HCC

p value was estimated by Kaplan-Meier method.

Table 2

Analysis of Association between Mitochondrial Haplogroups and HCC

mtDNA Haplogroup	Patients (n = 388)	Controls (n = 511)	OR (95% CI)	p Value
Macro-haplogroup
A	31 (8.0)	42 (8.2)	0.800 (0.465–1.375)	0.419
B	57 (14.7)	72 (14.1)	0.906 (0.584–1.406)	0.660
D	129 (33.2)	140 (27.4)	1.0
F	41 (10.6)	63 (12.3)	0.680 (0.422–1.095)	0.112
G	21 (5.4)	22 (4.3)	0.997 (0.511–1.947)	0.994
M7	24 (6.2)	35 (6.8)	0.713 (0.394–1.289)	0.263
M8	35 (9.0)	46 (9.0)	0.695 (0.414–1.165)	0.167
M9	3 (0.8)	6 (1.2)	0.503 (0.118–2.143)	0.352
M10	6 (1.5)	12 (2.3)	0.475 (0.169–1.339)	0.159
N9	9 (2.3)	35 (6.8)	0.241 (0.110–0.531)	0.000∗
Others	32 (8.2)	38 (7.4)	0.889 (0.514–1.540)	0.676
Age			1.001 (0.990–1.012)	0.876
Gender			0.304 (0.223–0.413)	0.000∗
Sub-haplogroup
B4	36 (9.3)	36 (7.0)	1.099 (0.626–1.928)	0.743
B5	16 (4.1)	28 (5.5)	0.776 (0.384–1.571)	0.482
CZ	25 (6.4)	37 (7.2)	0.637 (0.350–1.159)	0.140
D4	91 (23.5)	101 (19.8)	1.0
D5	36 (9.3)	38 (7.4)	1.091 (0.624–1.907)	0.759
F1	30 (7.7)	37 (7.2)	0.811 (0.455–1.445)	0.477
N9a	7 (1.8)	26 (5.1)	0.248 (0.101–0.610)	0.002∗
Y	2 (0.5)	9 (1.8)	0.253 (0.051–1.246)	0.091
Othersa	69 (17.8)	100 (19.6)	0.753 (0.488–1.162)	0.200
Age			1.001 (0.990–1.013)	0.795
Gender			0.307 (0.225–0.417)	0.000∗

p < 0.004 (0.05/14), adjusted p value with Bonferroni correction while 14 haplogroups were studied; values in parentheses are the percentage of samples. CI, confidence interval; OR, odds ratio.

Haplogroups with frequencies less than 5% in both patients and controls.

To ask whether HBV/HCV infection plays a role in haplogroup-associated HCC occurrence, we divided patients into HBV/HCV-positive and HBV/HCV-negative groups and compared them with the control group separately. As shown in Table S1, haplogroup N9 (OR 0.191, [95% CI 0.076–0.479], p < 0.001) and its sub-haplogroup N9a (OR 0.213, [95% CI 0.077–0.590], p = 0.003) are negatively associated with the occurrence of HBV/HCV-positive HCC. However, haplogroup N9 and its sub-haplogroup N9a are not associated with the occurrence of HBV/HCV-negative HCC (Table S2).

Taken together, we conclude that there is a significant negative correlation between haplogroup N9a and the prevalence of HBV/HCV-positive HCC in the study participants. Because only a small portion of patients (n = 43) with HCC were diagnosed without HBV/HCV infection, further studies with a greater number of HBV/HCV-negative HCC patients are needed to fully understand the role of mtDNA haplogroups in HBV/HCV-free HCC.

mtDNA SNP Analysis of Patients with HCC

To determine the prevalence of distinct mitochondrial SNPs (mtSNPs) in patient and control groups, we used ten SNPs that are preferentially expressed by distinct haplogroups (i.e., defining SNPs). To test whether the sample size we studied here was large enough, we determined a minimal sample size that provides sufficient statistical power to detect the OR = 0.5 or OR = 2 with the following conditions: average population minimum allele frequency (MAF) = 20% (average MAF in control group from this study is 31.24%) and case/control ratio = 1, p = 0.05, power = 90%. The calculated minimal sample size was 217 for OR < 0.5 and 180 for OR > 2, indicating that the sample size we used in this analysis (388) was qualified. After adjusting for age and sex, we found that the frequency of m.16257C > A, a defining mtSNP of the N9a haplogroup, was significantly lower in patients than in control participants (OR 0.290, [95% CI 0.123–0.685], p = 0.005) (Table 3). To exclude the false-positive probability of this significant association between m.16257C > A and HCC, a false-positive report probability (FPRP) analysis was used. As shown in Table S3, our results showed that with prior probability of 0.25 and 0.1, FPRP values were less than their significant predetermined value (0.2). This result suggested that our findings are deserving of attention. In addition, m.16189T > C, one of the common mtSNPs, was found higher in patients than in control participants (OR 1.497, [95% CI 1.115–2.010], p = 0.007); however, after applying the Bonferroni correction to adjust for ten simultaneous SNP statistical analysis, only p < 0.005 (0.05/10) was considered statistically significant (see Table 3). By dividing the patients into HBV/HCV-positive and HBV/HCV-negative groups, we found that both m.16189T > C (OR 1.486, [95% CI 1.137–1.943], p = 0.004) (Table S4) and m.16257C > A (OR 0.262, [95% CI 0.107–0.643], p = 0.003) (Table S4) are positively and negatively associated with the occurrence of HBV/HCV-positive HCC, respectively. However, m.16189T > C and m.16257C > A were found to not be associated with the occurrence of HBV/HCV-negative HCC (Table S5). The prevalence of other mtSNPs we evaluated, including the macro-haplogroup M-defining mtSNP, m.10400C > T, and the haplogroup D5a-defining mtSNP, m.10397A > G, were similar in patient and control groups regardless of the HBV/HCV infection status. However, further studies with a greater number of HBV/HCV-negative HCC patients are needed to fully understand the role of mtDNA in HBV/HCV-free HCC.

Table 3

Analysis of Association between Mitochondrial DNA SNPs and HCC

mtSNP	Patients (n = 388)	Controls (n = 511)	OR (95% CI)	p Value
m.150C > T	68 (17.5)	94 (18.4)	0.958 (0.671–1.367)	0.813
m.249delA	50 (12.9)	85 (16.6)	0.703 (0.477–1.037)	0.076
m.10397A > G	38 (9.8)	38 (7.4)	1.426 (0.875–2.324)	0.155
m.10398A > G	254 (65.5)	316 (61.8)	1.189 (0.895–1.580)	0.232
m.10400C > T	229 (59.0)	278 (54.4)	1.154 (0.875–1.520)	0.310
m.16129G > A	63 (16.2)	106 (20.7)	0.702 (0.492–1.001)	0.051
m.16189T > C	138 (35.6)	147 (28.8)	1.497 (1.115–2.010)	0.007
m.16223C > T	238 (61.3)	293 (57.3)	1.066 (0.806–1.410)	0.652
m.16257C > A	7 (1.8)	26 (5.1)	0.290 (0.123–0.685)	0.005∗
m.16362T > C	188 (48.5)	214 (41.9)	1.293 (0.983–1.702)	0.066

p < 0.005 (0.05/10), adjusted p value with Bonferroni correction while 10 SNPs were studied. CI, confidence interval; OR, odds ratio. Values in parentheses are the percentage of samples.

These data indicate that m.16189T > C and m.16257C > A are positively and negatively associated with the prevalence of HBV/HCV-positive HCC, respectively. The associations between mtSNPs and HBV/HCV-negative HCC are currently not known, and further studies with greater number of HBV/HCV-negative HCC patients are needed.

Evaluation of Clinical Characteristics Associated with Distinct mtDNA Haplogroups

To evaluate whether distinct mtDNA haplogroups were associated with specific clinical characteristics of HCC, we analyzed several clinical characteristics including alpha-fetoprotein (AFP) levels, tumor size, and biomarkers of liver function such as aspartate aminotransferase (AST), alanine aminotransferase (ALT), serum albumin (ALB), and total bilirubin (TBIL) levels. We did not correct for differences in HBV infection status and tumor metastasis among haplogroups, because only a small percentage of patients tested negative for HBV infection (15.49%) and tumor metastasis (7.94%). With our original number of recruited participants, dividing haplogroups further based on these two characteristics would diminish the power of our study. A recent study found that liver function biomarkers were elevated in patients with HCC, suggesting that biomarkers of liver function may be closely associated with the risk of developing HCC. However, using multivariate logistical regression analysis with adjustment for age, sex, and HBV/HCV infection, we found no correlation between mtDNA haplogroups and AST (Table S6), ALT (Table S7), ALB (Table S8), or TBIL (Table S9). Next, we analyzed the potential association between mtDNA haplogroups and AFP expression in patients with HCC. Based on the criterion that AFP > 400 ng/mL is a diagnostic index of HCC, patients were divided into high AFP expression group (>400 ng/mL) and low AFP expression group (<400 ng/mL); however, haplogroups including N9a were not associated with the level of AFP expression level (Table 4). By dividing the patients into group with large tumor size (diameter > 5 cm) and group with small tumor size (diameter > 5 cm), no haplogroups were found associated with tumor size (Table 5). Notably, we found that mean AFP levels and tumor size were lowest in patients with haplogroup N9a and were top ranked in patients with haplogroup D5 and B5, respectively (Figures S1A and S1B). Taken together, our results suggest that mtDNA haplogroup N9a is not associated with clinical characteristics of HCC.

Table 4

Association between mtDNA Haplogroup and Serum AFP Level

mtDNA Haplogroup	AFPa		Multivariate
	High (n = 117) [n (%)]	Low (n = 256) [n (%)]	OR (95% CI)	p Value
A	11 (9.4)	19 (7.4)	1.005 (0.417–2.418)	0.992
B4	16 (13.7)	18 (7.0)	1.651 (0.728–3.743)	0.230
B5	3 (2.6)	11 (4.3)	0.493 (0.124–1.965)	0.316
CZ	5 (4.3)	20 (7.8)	0.418 (0.140–1.248)	0.118
D4	31 (26.5)	59 (23.0)	1.0
D5	10 (8.5)	24 (9.4)	0.743 (0.309–1.790)	0.508
F1	12 (10.3)	15 (5.9)	1.536 (0.630–3.743)	0.345
G	5 (4.3)	15 (5.9)	0.729 (0.238–2.230)	0.579
M7	4 (3.4)	20 (7.8)	0.356 (0.110–1.157)	0.086
N9a	1 (0.9)	6 (2.3)	0.354 (0.040–3.128)	0.350
Y	0 (0.0)	1 (0.4)	0.000 (0.000)	1.000
Others	19 (16.2)	48 (18.8)	0.787 (0.391–1.584)	0.502
Age			0.966 (0.945–0.987)	0.002∗
HBV/HCV			1.179 (0.542–2.565)	0.679
Gender			0.826 (0.452–1.510)	0.534

a15 patients were excluded and AFP level and HBV/HCV infection status were not available for 14 patients and 1 patient, respectively.

p < 0.004 (0.05/12), adjusted p value with Bonferroni correction while 12 haplogroups were studied; CI, confidence interval; OR, odds ratio.

Table 5

Association between mtDNA Haplogroup and Tumor Size

mtDNA Haplogroup	Tumor Sizea		Multivariate
	Large (n = 142) [n (%)]	Small (n = 170) [n (%)]	OR (95% CI)	p Value
A	15 (10.6)	12 (7.1)	1.470 (0.602–3.587)	0.397
B4	9 (6.3)	21 (12.4)	0.475 (0.191–1.184)	0.110
B5	9 (6.3)	3 (1.8)	3.599 (0.884–14.653)	0.074
CZ	7 (4.9)	12 (7.1)	0.626 (0.221–1.773)	0.378
D4	37 (26.1)	42 (24.7)	1.0
D5	11 (7.7)	17 (10.0)	0.761 (0.311–1.860)	0.549
F1	9 (6.3)	12 (7.1)	0.806 (0.302–2.150)	0.667
G	6 (4.2)	7 (4.1)	1.059 (0.323–3.473)	0.925
M7	7 (4.9)	13 (7.6)	0.614 (0.218–1.724)	0.354
N9a	3 (2.1)	4 (2.4)	0.819 (0.167–4.007)	0.805
Y	0 (0.0)	1 (0.6)	0.000 (0.000)	1.000
Others	29 (20.4)	26 (15.3)	1.302 (0.647–2.622)	0.460
Age			0.981 (0.960–1.003)	0.089
HBV/HCV			0.687 (0.330–1.428)	0.314
Gender			0.720 (0.402–1.289)	0.269

CI, confidence interval; OR, odds ratio.

a76 patients were excluded and tumor size and HBV/HCV infection status were not available for 72 patients and 4 patients, respectively.

N9a Cybrids and Non-N9a (B5, D4, and D5) Cybrids Feature Distinct Transcriptome Profiles and Neoplastic Activity

To investigate the contribution of N9a in tumorigenesis, we generated five cybrid cells by fusion of mtDNA-lacking Rho zero human osteosarcoma 143B cells with platelet containing mtDNA haplogroup N9a (N9a10a and N9a1) and non-N9a cybrids (B5, D4, and D5), respectively. In vitro neoplastic assay revealed that growth rate of N9a cybrids was significantly slower than non-N9a cybrids (Figure 3A). Moreover, the number of colonies formed in N9a cybrids was significantly less than non-N9a cybrids (Figure 3B). These data suggest that cells with haplogroup N9a have lower tumorigenic activity than non-N9a cybrids.

Figure 3

Haplogroup N9a Have Lower Tumorigenic Activity than Non-N9a (B5a1, D5a2, D4g2) Cybrids

(A) Cell proliferation curve of N9a and non-N9a cybrids (n = 3). (B) Colony formation assay of N9a and non-N9a cybrids (n = 3). Data are presented by means ± SEM, and ***p < 0.001.

Next, we sought to identify the molecular mechanisms responsible for the negative correlation between haplogroup N9a and HCC. We analyzed the transcriptome of two N9a cybrids (N9a10a and N9a1) and three non-N9a cybrids (B5, D4, and D5). Multiple genes were differentially expressed between N9a cybrids and non-N9a cybrids; the transcription of 2,748 genes increased and the transcription of 3,185 genes decreased (p < 0.05, |log2FC| > 1) (Figure 4A). Using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database as a reference, gene set enrichment analysis (GSEA) revealed that 49 pathways were differentially regulated in N9a cybrids compared to non-N9a cybrids (Data S1). Nearly 50% of the KEGG pathways identified by GSEA were associated with cellular metabolism (Figure 4B), and most of these, including fatty acid metabolism and fatty acid synthesis, were downregulated. Further, we found that increased activity in the oxidative phosphorylation (OXPHOS) pathway was positively correlated with the prevalence of haplogroup N9a in study participants (Figure 4B). However, in a previous study, we demonstrated that increased expression of nuclear DNA-encoded but not mtDNA-encoded OXPHOS subunits in N9a cybrids may compensate for diminished mitochondrial function (Figures 4C and 4D), and suggested that the N9a haplogroup exhibited lower mitochondrial OXPHOS activity than non-N9a cybrids. We performed GSEA using the Reactome and CGP (chemical and genetic perturbations) gene sets from the Molecular Signatures Database (MSigDB) collection. Moreover, positive regulatory machinery of cell-cycle progressions was less activated in N9a cybrids than non-N9a cybrids. Alternatively, cell-cycle checkpoint pathways (including both p53-dependent and p53-independent checkpoints) were more activated in N9a cybrids than non-N9a cybrids (Figure 4E). Overall, the analysis found that the genetic profiles of non-N9a cybrids positively correlated with liver cancer invasion, progression, and recurrence, whereas the genetic profiles of N9a cybrids negatively correlated with liver cancer (Figure 4F).

Figure 4

N9a and Non-N9a (B5a1, D5a2, D4g2) Cybrids Feature Distinct Transcriptome Profiles

(A) Heatmap of RNA-seq data showing the differentially expressed genes (DEGs) between N9a and non-N9a cybrids (p < 0.05, |log2FC| > 1). (B) GSEA enrichment using KEGG database showing metabolic pathways are differently activated between N9a and non-N9a cybrids. (C and D) Heatmap showing transcriptional changes of nuclear-encoded (C) and mitochondrially encoded (D) OXPHOS subunits in N9a and non-N9a cybrids. (E) GSEA enrichment of cell-cycle related events in N9a and non-N9a cybrids using KEGG, Reactome, and GCP from MsigDB. (F) GSEA enrichment of liver cancer related events in N9a and non-N9a cybrids using GCP.

Discussion

The association of mtDNA haplogroups with HCC in China varied by geographical region. In northern China, haplogroup M and its defining SNP, m.T489C, were positively associated (albeit with marginal significance) with later age-at-onset of HCC, although haplogroup M7 was negatively associated with HCC in eastern China. In this study, we found that haplogroup N9a was negatively associated with HCC in northern China, but neither haplogroup M7 nor macro-mtDNA haplogroup M (defined by m.10400C > T) was associated with HCC in the same region.

Previous studies found that haplogroup N9a was negatively associated with type 2 diabetes (T2D) and other metabolic syndromes in Japan and Korea; however, this association was not found in Taiwan. Recently, we identified N9a as a risk factor for T2D in the southern Chinese population, which suggests that haplogroup N9a may play a geographic specific role in the development of metabolic disease. While HCC is positively co-related with metabolic diseases including non-alcoholic fatty liver disease (NAFLD) and T2D,29, 30 it is reasonable that the haplogroup N9a may play the same role in the incidence of HCC. Because other factors may differentially affect the role of haplogroup N9a on T2D in different regions,11, 27, 28 it is difficult to speculate whether haplogroup N9a plays a consensus or opposite role in HCC and T2D in northern China. Further studies on the association between haplogroup N9a and metabolic diseases in northern China may help. Notably, we found that mean AFP and mean tumor volume are top ranked in patients with haplogroup D5 and haplogroup B5, indicating that haplogroup B5 and D5 may play a role in the disease progression of HCC. Interestingly, our previous report showed that cybrid with haplogroup D5 is positively associated with cancer progression, such as cancer cell reproductive viability, transforming potential, and migration.10, 31 Although associations between haplogroup B5 and cancers were not previously reported, a positive association between haplogroup B5 and Alzheimer’s disease susceptibility was reported. These suggested that haplogroup B5 and D5 may play a role in the progression of HCC. However, due to the limit of sample size and lack of follow-up of patients with HCC, further studies with large sample size and clinical characteristics with time trends are needed to fully address the association between mtDNA haplogroups and clinical features.

Both increased and decreased mitochondrial OXPHOS functions were shown to promote tumorigenesis in different cancers with distinct mechanisms.32, 33 For some mtDNA haplogroups (e.g., haplogroup D5 in breast cancer), lower mitochondrial function may promote tumorigenesis due to altered mitochondrial signals including reactive oxygen species and NAD+/NADH, which are responsible for overactivation of mitochondrial retrograde signaling. Although haplogroup N9a was positively associated with OXPHOS function in GSEA using microarray data in one study, this conclusion was not supported by the following mitochondrial function study. In our recent study, we confirmed that N9a is positively associated with OXPHOS-related genes expression; however, we found that mitochondrial function in N9a cybrids was lower than in non-N9a cybrids, suggesting that the positive correlation between OXPHOS GSEA and haplogroup N9a may be a means of compensating for reduced mitochondrial function in N9a cybrids. In this study, a positive association between OXPHOS GSEA and haplogroup N9a was confirmed further by using two N9a cybrids and three non-N9a cybrids (B5, D4, and D5) (Figure 4B). Therefore, it is likely that haplogroup N9a plays a protective role in HCC development by decreasing mitochondrial function. Increased mitochondrial function was shown to promote tumorigenesis by generating ATP for cancer-related kinase activation, supplying tricarboxylic-acid-cycle-derived amino acid for endothelial cell proliferation during angiogenesis, and driving epithelial-mesenchymal transition progress for metastasis. In this study, we did see an increased metastatic activity in non-N9a cybrids compared with N9a cybrids by using MSigDB-based GSEA (Data S2). Furthermore, GSEA using KEGG revealed that half (24/49) of enriched pathways are metabolic pathways (Data S1). Given that change of metabolic profiles is one of the key features in tumorigenesis, we believe that N9a and non-N9a haplogroups play different role in HCC incidence. However, we do not know how cellular metabolism pathways are regulated by different performance of mitochondrial function in N9a and non-N9a haplogroup containing cells. Further studies on the causal relationship between OXPHOS function and cellular metabolism pathways are necessary to completely reveal the protective role of haplogroup N9a in HCC.

In conclusion, we identified a negative correlation between mtDNA haplogroup N9a and HCC prevalence and tumor volume in patients with HCC. We also demonstrated that cells containing haplogroup N9a and non-haplogroup N9a exhibit distinct metabolic profiles, which may contribute to the negative correlation between N9a and the probability of developing HCC. Additional studies are needed to fully characterize the metabolic basis for the effects of haplogroup N9a on HCC incidence and progression.

Materials and Methods

Study Participants

From March 2018 to July 2019, we recruited 388 patients with HCC (mean ± SD 55.68 ± 0.55, median age 55, range 12–88) from the Henan Provincial People’s Hospital (n = 209, mean ± SD 56.18 ± 0.79, median age 57, range 12–88), the First Affiliated Hospital of Zhengzhou University (n = 115, mean ± SD 55.23 ± 0.94, median = 55, range 27–77), and the Henan Provincial Cancer Hospital (n = 64, mean ± SD 54.86 ± 1.32, median = 53.5, range 34–73) to participate in the study. Each diagnosis of HCC was confirmed histologically, either immediately after surgery or with a subsequent liver biopsy. Meanwhile, we recruited 511 geographically matched, asymptomatic control participants (mean ± SD 56.19 ± 0.59, median = 58, range 19–87) with no history of cancer or of other diseases associated with mitochondrial defects (e.g., diabetes and hypertension), from The Physical Examination Center of Henan Provincial People’s Hospital from May 2019 to July 2019. In the patient group, 338 patients are HBV or HCV positive (320 patients are HBV positive, 16 patients are HCV positive, and 2 patients are HBV+HCV positive), 43 patients are HBV/HCV negative, and HBV/HCV infection history was missing in 7 patients. All control subjects are HBV/HCV negative. A detailed characteristic of age and sex of the patients and control subjects was described in Table S10. Informed consent was obtained from the subjects enrolled in this study, and the study was approved by the Ethics Committee of Wenzhou Medical University and all three hospitals we mentioned above.

mtDNA Sequencing and Genotyping

Genomic DNA was extracted from peripheral blood using a standard SDS lysis protocol. Two primer pairs were designed to amplify mtDNA fragments, which contain major diagnostic SNPs of the Asian mtDNA haplogroup: L15975F: 5′-CTCCACCATTAGCACCCAAAGC-3′ and H794R: 5′-AGGCTAAGCGTTTTGAGCTG-3′; L9967F: 5′-TCTCCATCTATTGATGAGGGTCT-3′ and H10858R: 5′-AATTAGGCTGTGGGTGGTTG-3′. An additional sequencing primer, 299F 5′-GGTGGAAATTTTTTGTTATG-3′, was designed to detect a potential poly C gap between mt16184 and mt16193. PCR was performed using a PCR Amplification Kit (TaKaRa, Tokyo, Japan) with a S1000 Thermal Cycler (Bio-Rad, Hercules, CA, USA) under the following conditions: pre-denaturation at 95°C for 5 minutes, followed by 35 cycles of (94°C for 30 s, 57°C for 30 s, and 72°C for 40 s), with a final extension at 72°C for 4 min. Sanger sequencing was performed using ABI 3730XL (Thermo Fisher Scientific, Waltham, MA, USA). SNPs of each participant were identified by comparing sequences with the revised Cambridge Reference Sequence (rCRS) using CodonCode Aligner 3.0.1 (CodonCode Corporation, Centerville, VA, USA). A mtDNA haplogroup was identified for each participant using MitoTool, and visually confirmed by comparing SNPs from the D-loop of ND3 and ND4L, with the diagnostic SNPs from the east Asian mtDNA haplogroup tree.37, 38, 39 If haplogroup identification was still inconclusive, additional information was obtained by restriction fragment-length polymorphism (RFLP) analysis at sites: 3,010 (BccI), 4,833 (HhaI), 5,178 (AluI), and 9,824 (HinfI) and by testing for the 9 bp deletion at the COII-tRNAlys junction. mtDNA sequence variations and RFLP results for each individual can be found in Data S3.

In Vitro Tumorigenesis Assay

For cell proliferation assay, a total of 100,000 cells were cultured in 6-well plate, and the number of cells was counted for 24 h, 48 h, 72 h, and 96 h. For colony formation assay, 1,000 cells were cultured in 6-well plate for 2 weeks and then fixed with 4% paraformaldehyde (30525-89-4,Yonghua, Shanghai, China) for 30 min, followed by staining with 0.1% crystal violet solution (C0121, Beyotime, Shanghai, China) for 30 min, and then washed three times with PBS buffer. ImageJ v2.4.1.7 was used to analyze the number of clones.

RNA Sequencing and Analysis of Gene-Expression Data

Two N9a cybrids (N9a10a and N9a1) and three non-N9a (B5, D4, and D5) cybrids were generated in our previous study. To evaluate the differences between the gene-expression profiles of N9a cybrids and non-N9a cybrids, we isolated total RNA from each group of cybrids (in triplicate) using a RNeasy Mini Extraction kit (QIAGEN, Valencia, CA, USA), and mRNA from 20 μg of the total RNA was purified using poly-T-attached magnetic beads. After fragmenting the mRNA, first-strand cDNA was synthesized and sequenced using an Illumina HiSeq 2000 platform (Illumina, San Diego, CA, USA) as described previously. Clean reads were obtained by deleting adaptor-only sequences and low-quality sequences. Sequence comparison was carried out using STAR (v2.5.1b) against 1000 Genomes Build 37 Decoy 5. The number of annotated clean reads of each gene was analyzed and normalized as reads per kilobase per million reads (RPKM). Genes were considered differentially expressed if their transcript level changed at least 2-fold (p < 0.05) compared with controls. GSEA was performed using the KEGG and Reactome databases as references. To evaluate the enrichment of liver cancer-related pathways, we performed GSEA using gene sets of CGP from the MSigDB as references. The minimum and maximum genes number were 30–200 for GSEA. Pathways were considered significantly enriched with false discovery rate < 0.25 and p < 0.05.

Statistical Analysis

In this case-control study, statistical analysis of mtDNA haplogroup and HCC was performed in three levels: macro-haplogroups M and N, major branches of haplogroups M and N according to the phylogenetic tree of Han Chinese, and haplogroup with frequency of >5% in either the control group or the patient group; these were included to evaluate the potential association between common mtDNA haplogroups and HCC. All other haplogroups with frequencies of <5% were collectively evaluated as “other haplogroups.” PCA of the mtDNA haplogroup frequencies was conducted by the R package to show the clustering pattern of patients with HCC and control subjects, as well as reported Han Chinese population across China (Table S11). The minimal sample size study was determined by using Quanto software in the study of associations between mtDNA SNPs and patients with HCC. FPRP analysis was calculated by the R package. Multivariate logistical regression analysis was applied to adjust for disparities in age and gender between the patient group and control group. We found that for each mtDNA haplogroup, HCC was a dependent variable, whereas age, sex, and genotype were independent variables. Because there was high numerical variability among the 14 mtDNA haplogroups included in the analysis (i.e., A, G, M7, M9, M10, B4, B5, CZ, D4, D5, F1, N9a, Y, and “other haplogroups”), we applied dummy coding. A Bonferroni correction indicating p < 0.004 (0.05/14) was considered statistically significant when evaluating these haplogroups. To examine the relationships among liver function biomarkers and the mtDNA haplogroups, we applied multivariate logistical regression analysis to adjust for risks associated with age, sex, and HBV infection status. One-way ANOVA was used to analyze the differences in mean AFP expression and volume among ten mtDNA haplogroups. In vitro tumorigenesis assays were analyzed using two-tailed Student’s t test. All statistical analyses were performed using SPSS 21.0 (IBM, Armonk, NY, USA) unless otherwise indicated.

Author Contributions

L.S., M.J., and H.F. conceived and designed the study. S.H., M.L., J.W., Y.Z., and J.L. collected the samples, performed the PCR, and aligned sequences. Q.Z. and J.L. performed the RNA-sequencing experiment and the following bioinformatical analysis. L.S., M.J., and H.F. supervised the work and wrote the manuscript.