Genetic Alterations in Mitochondrial DNA Are Complementary to Nuclear DNA Mutations in Pheochromocytomas

Simple Summary Mitochondrial DNA (mtDNA) alterations have been reported to play important roles in cancer development and metastasis. However, there is scarce information about pheochromocytomas and paragangliomas (PCCs/PGLs) formation. To determine the potential roles of mtDNA alterations in PCCs/PGLs, we analyzed a panel of 26 nuclear susceptibility genes and the entire mtDNA sequence of 77 human tumors, using NGS. We also performed an analysis of copy-number alterations, large mtDNA deletion, and gene/protein expression. Our results revealed that 53.2% of the tumors harbor a mutation in the susceptibility genes and 16.9% harbor complementary mitochondrial mutations. Large deletions and depletion of mtDNA were found in 26% and 87% of tumors, respectively, accompanied by a reduced expression of the mitochondrial biogenesis markers (PCG1α, NRF1, and TFAM). Furthermore, P62 and LC3a gene expression suggested increased mitophagy, which is linked to mitochondrial dysfunction. These finding suggest a complementarity and a potential contributing role in PCCs/PGLs tumorigenesis. Abstract Background: Somatic mutations, copy-number variations, and genome instability of mitochondrial DNA (mtDNA) have been reported in different types of cancers and are suggested to play important roles in cancer development and metastasis. However, there is scarce information about pheochromocytomas and paragangliomas (PCCs/PGLs) formation. Material: To determine the potential roles of mtDNA alterations in sporadic PCCs/PGLs, we analyzed a panel of 26 nuclear susceptibility genes and the entire mtDNA sequence of seventy-seven human tumors, using next-generation sequencing, and compared the results with normal adrenal medulla tissues. We also performed an analysis of copy-number alterations, large mtDNA deletion, and gene and protein expression. Results: Our results revealed that 53.2% of the tumors harbor a mutation in at least one of the targeted susceptibility genes, and 16.9% harbor complementary mitochondrial mutations. More than 50% of the mitochondrial mutations were novel and predicted pathogenic, affecting mitochondrial oxidative phosphorylation. Large deletions were found in 26% of tumors, and depletion of mtDNA occurred in more than 87% of PCCs/PGLs. The reduction of the mitochondrial number was accompanied by a reduced expression of the regulators that promote mitochondrial biogenesis (PCG1α, NRF1, and TFAM). Further, P62 and LC3a gene expression suggested increased mitophagy, which is linked to mitochondrial dysfunction. Conclusion: The pathogenic role of these finding remains to be shown, but we suggest a complementarity and a potential contributing role in PCCs/PGLs tumorigenesis.


Introduction
Pheochromocytomas (PCCs) and paragangliomas (PGLs) are rare endocrine tumors that arise from the adrenal medulla and in paraganglia of the autonomous nervous system, respectively. Up to 40% of pheochromocytomas are hereditary and caused by germline mutations in well-known cancer susceptibility genes (SDHx, VHL, EGLN1/PHD2, EPAS1/HIF2A, KIF1Bβ, MAX, MEN1, NF1, RET, TMEM127, RAS, NF1, FGFR1, ATRX, etc.) which are involved in interconnecting pathways [1,2]. Sporadic pheochromocytomas can have somatic mutations in one of these hereditary genes, indicating that they may be the main drivers. Together, germline and somatic mutations in susceptibility genes are found in approximately 60% of all PCCs/PGLs [3,4], but the cause of many sporadic PCCs/PGLs is still unknown [5]. With regard to gene-expression patterns, at least two main signatures are distinguished, clusters 1 and 2. Cluster 1 tumors (having mutations in, for example, VHL and SDHA/B/C/D/AF2) are characterized by increased expression of genes involved in (pseudo) hypoxia, mitochondrial electron transport chain, Krebs cycle, cell proliferation, and angiogenesis. Cluster 2 tumors (having mutations in, for example, RET, NF1, MAX, and TMEM127) are characterized by an increased expression of genes involved in cell proliferation, protein synthesis, and kinase signaling. Recently, a set of 46 genes was used by Flynn et al. (2016) to develop a Pheo-type gene-expression profiling, where the samples were separated into two main subtypes: pseudohypoxia and RTK/Ras driven tumors. Pseudohypoxia clustering is defined by overexpression of the transcription factor, TRIB3, and the cell adhesion gene, DSP, for SDHx tumors and overexpression of angiogenesis-associated genes (i.e., ETS1, CD34, and AQP1) for VHL tumors. RTK clustering is based on the overexpression of a chromaffin cell differentiation signature (i.e., PNMT and RET) within NF1, RET, HRAS, and TMEM127 tumors and underexpression of genes associated with the MAP kinase pathway (i.e., SPRY4 and DUSP6) and metabolism (i.e., CSGALNACT1 and PFKFB3) for MAX-like tumors [6].
The nuclear genes representing cluster 1 have a main role in the pathogenesis of PCCs/PGLs and other types of cancer by affecting the cellular metabolism and mitochondrial function. However, cancer pathogenesis commonly involves the accumulation of several genetic alterations, not only in the nuclear genome but also in the mitochondrial genomes. In the past few years, several studies have reported a significant number of mtDNA mutations, including point mutations, deletions, and insertions, in a variety of human malignancies, including head and neck [7], lung [8], and breast [9]. It is believed that mtDNA mutations and copy-number alterations are associated with tumor malignancies, because of the high levels of reactive oxygen species (ROS) produced during oxidative phosphorylation (OXPHOS) and oxidative mtDNA damage, along with a less efficient mtDNA repair system [10,11]. The depletion of mtDNA has also been identified in various cancers [12] and has been suggested as a potential prognostic marker in cancers of the bladder, breast, kidney, head and neck, esophagus, and liver [13,14].
Recent publications have shown that catecholamine metabolism is fundamental to mtDNA integrity and mitochondrial function [15,16], but mitochondrial genome analysis has not been investigated in PCCs/PGLs yet.
In this study, we investigated mutational status/mutations in a panel of 26 susceptibility genes and the mitochondrial genome in PCCs/PGLs tumors by performing an extensive analysis of mutations, using next-generation sequencing, mitochondrial copy-number alterations, large deletions, and gene and protein expression.

Patients and Tumors
This study includes 74 pheochromocytomas and 3 paragangliomas from 77 patients having a sporadic disease, i.e., without family history and any syndromic features (see Supplementary Table S1). Fifty-four tumors were obtained from the Brabois Hospital, Nancy, France, and the remaining 23 tumors were obtained from Linkoping University Hospital, Sweden. All tumors were fresh-frozen. As a control, peripheral blood was obtained from 22 Swedish patients, and normal fresh-frozen tissues were obtained from 50 French patients. All normal tissues from the adrenal medulla were carefully dissected from the tumors during the surgery, at a distance normally larger than 1 mm, followed by tumor dissection and macroscopic aspect. All samples were histologically confirmed as PCCs/PGLs using World Health Organization (WHO) criteria and the distinction between malignant and benign PCCs/PGLs was evaluated according to the Endocrine Society Clinical Guidelines Subcommittee (CGS) criteria [17]. All tissue samples were handled following a Standard Operation Procedure. In short, resection specimens were stored in labeled cryovials and snap-frozen in liquid nitrogen. The time laps between sample resection and freezing were less than 15 min. All samples were collected and studied (Appendix A) with informed consent and approval from the local ethics committees (Dnr 2010/40-31, Dnr 2015/175-32 Linköping University, Sweden; DR-2016-346, CHRU de Nancy, France).

Mitochondrial Copy-Number Variation and DNA Large Deletion
Mitochondrial DNA-copy number analysis was performed by digital droplet PCR (ddPCR), as described previously by Memon et al. [33], using probes targeting the mitochondrially encoded NADH dehydrogenase 1 (MT-ND1) (assay ID: dHsaCPE5029120, Bio-Rad) gene labeled with FAM fluorophore (mitochondrial DNA) and the eukaryotic translation initiation factor 2C, 1 (EIF2C1) (assay ID: dHsaCP1000002, Bio-Rad), and ribonuclease P/MRP 30 kDa subunit (RPP30) (assay ID: dHsaCP1000485, Bio-Rad) genes (nuclear DNA) labeled with HEX. Normal adrenal medulla tissues from the French cohort were used as a control group to estimate the mitochondrial copy-number variation in both tumor cohorts.
Long-range PCR was performed to detect mtDNA deletions, using Long PCR Enzyme Mix, according to manufacturer's instructions (LA Taq DNA polymerase # RR002A, Takara). To amplify the major arc and minor arc, we used two pair of primers, F: 5 TGGCTCCTTTAACCTCTCCA 3 ; R: 5 AGGCTAAGCGTTTTGAGCTG 3 and F: 5 CTCCT-CAAAGCAATACACTG 3 ; and R: 5 AAGGATTATGGATGCGGTTG 3 , respectively. The amplified fragments of the 77 samples and their corresponding controls were analyzed with the fragment analyzer (5200 Fragment Analyzer System, Agilent, Santa Clara, CA, USA).

Sequencing of Nuclear Genes Involved in Mitochondrial DNA Integrity
The coding regions of POLG1 (NM_002693), C10ORF2 (Twinkle) (NM_021830), and TFAM (NM_003201) genes were sequenced with the dideoxy nucleotide termination method (Sanger sequencing). PCR was performed on cDNA, using primers that amplify the complete coding region of the transcript. All identified cDNA variants were confirmed on the corresponding genomic DNA. For the analysis of DGUOK (NM_080916) and MPV17 (NM_002437), all exons were amplified with intronic primers and sequenced with the Sanger methodology on an ABI 3500 sequencer (Applied Biosystems) (Appendix C) (see Supplementary Table S3 for primers).

Microarray Gene Expression Analysis
Samples used for microarray analysis had RNA integrity number (RIN) values ≥ 8.5. Thirty-nine samples were excluded from this analysis, due to the lack of samples or to low RIN values. Generated sense-strand cDNA, using the Ambion WT Expression Kit (Life Technologies, Carlsbad, CA, USA), was fragmented, labeled, and hybridized to Human Gene 1.0 ST arrays (Affymetrix, Santa Clara, CA, USA), according to the manufacturer's protocols. The arrays were washed, stained, and scanned in a GeneChip Scanner 3000 7G (Affymetrix, Santa Clara, CA, USA). CEL files were analyzed with Transcriptome Analysis Console (TAC) Software V4.0 (Affymetrix). Genes with p-values < 0.05, FDR-corrected p-values < 0.25, and fold change <−2 or >2 were considered significantly regulated genes.
Since we were analyzing tumors from two different populations, we are aware of the possibility of different handling issues, even though the sample handling was performed carefully and similarly. The gene-expression patterns were analyzed separately by performing Gene Set Enrichment Analysis (GSEA), using GSEA software V4.0.3 (Broad institute) (accessed on 27 September 2021) and referring to the Molecular Signatures database (MSigDB) for the "C2: Canonical pathways" gene sets.
Given their involvement in mitochondrial biogenesis and autophagy, P62, LC3, PCG1α, and NRF1 were explored. GUSB was used as an endogenous control gene. RT-PCR analysis was performed by using SYBR Green, with which we used 10 ng of cDNA; all reactions were performed in duplicates, using the 7500 fast real-time PCR system (Applied Biosystem, Foster City, CA, USA) (see Supplementary Table S3 for primers).

Mass Spectrometry and Protein Analysis
The protein samples (mitochondria) were reduced by 5 mM of ditiothreitol, alkylated with 10 mM iodoacetamide, and incubated at 37 • C overnight with mass-spectrometrygrade trypsin (Thermo-Fisher, USA). The obtained peptides were desalted on Pierce™ C18 Tips (Thermo-Fischer, USA), according to manufacturer's instruction. The desalted peptides were vacuum-dried, diluted with 12 µL 0.1% formic acid, and used for mass spectrometry analysis (see Appendix D for further information). Protein-enrichment analysis to identify related pathways and their functional connectivity was performed by using ProtExA (http://bioinformatics.cing.ac.cy/protexa/) (accessed on 27 October 2021), a web tool for the statistical and functional analysis of protein-expression datasets [37].

Statistical Analysis
Generation of Receiver Operating characteristic (ROC) curve in SPSS statistical package (Version 25.0), with the area under curve (AUC), cutoff, sensitivity, and 1-specificity values, served to set mtDNA content cutoff value, discriminating between normal and tumor tissue (Appendix D).
For genes and proteins expression, the quantitative differences were statistically analyzed by using Student's t-test; differences with p-values lower than 0.05 were considered statistically significant. Graphs were generated by using GraphPad Prism software V5.0.

Novel Mutations in PCCs/PGLs Susceptibility Genes
To examine mutations in these two cohorts, we designed a gene panel comprising all known and PCCs/PGLs-implicated susceptibility genes. In this design, 26 genes were covered by 706 probes. Of these, 99.8% of the enriched regions yielded sequence reads with a mean depth of 1000× per target and sample. All 41 identified mutations were confirmed with Sanger sequencing, and their status was checked in corresponding normal tissues (blood for all Swedish samples or non-tumor adrenal medullas for all French samples). Our analysis of paired normal DNA revealed that 25 mutations were somatic (61%), 10 were germline (24.4%), and six could not be assessed due to the lack of normal tissue (14.6%). Bioinformatic analysis suggested that 41 tumors presented mutations which are most likely the drivers, with seven mutations being novel. Three mutations in the NF1; a deletion of 12-bp in exon 27 and two frameshift mutations, c.5704_5705insC (L1902fs) and c.1904_1907 del (P635fs); two mutations in EGLN1, a frameshift mutation c.607_619 del(N203fs), and a non-sense mutation c.153G>A (W51X). Two tumors displayed frameshift and nonsense mutations in SDHD (c.205_217 del (E69fs)) and CSDE1 (c.1660C>T (R554X)), respectively. In total, 53.2% of all PCCs/PGLs had mutations in the targeted susceptibility genes, leaving many tumors without a genetic explanation (Table 1) ( Figure 1A). Table 1. Pathogenic mutations identified among 41 patients with pheochromocytomas and paragangliomas from Nancy, France, and from Linköping, Sweden.

Sample Gene Mutation Protein Variation Mutation Status Previously Reported
French cohort      After classifying the samples according to mRNA expression subtypes, we added gene-expression data for the mtOXPHOS genes ( Figure 1). Most of the mtOXPHOS genes are generally upregulated over all Pheo-Type subgroups, except for MT-ND2/ND6 genes, which are downregulated, indicating probably a mitochondrial impairment in complex I among all PCCs/PGLs.

Molecular Profiling and Gene Expression Patterns
Our analysis of gene expression patterns revealed a difference and separated the Swedish and French cohorts into two groups. To better infer the biological differences, the two cohorts were analyzed separately, and a gene-set enrichment analysis was performed. Our analysis of gene expression patterns revealed a difference and separated the Swedish and French cohorts into two groups. To better infer the biological differences, the two cohorts were analyzed separately, and a gene-set enrichment analysis was performed. GSEA determined the functional categories of significantly affected genes and yielded 14 gene sets that are enriched in the Swedish cohort: 12 gene sets are downregulated and associated mainly with genes involved in protein metabolism and RNA translation, and two upregulated gene sets associated with GPCR (G protein coupled receptors) signaling (Supplementary Table S4 and Figure S1).

Novel Mitochondrial Mutations
Massive parallel sequencing of mtDNA revealed several variations in both known and novel polymorphic sites and allowed us to classify the 77 PCCs/PGLs patients into different haplogroups with different distributions separating the French and Swedish cohorts, but no significant association was found between mtDNA haplogroup bins and PCCs/PGLs (Supplementary Table S5 and Figure S2).
The POTTER software was used to predict the effect of the novel two nonsense and frameshift somatic mutations on the secondary structures of the ND1, ND6 (complex I, NADH Coenzyme Q oxidoreductase), and COXIII (complex IV, Cytochrome c Oxidase) proteins, respectively (Supplementary Figure S3A). The screening of genetic variants in mitochondrial genes revealed eight novel missense mutations, including three germline mutations (m.8466A>T, m.12068A>G, and m.15471T>C) in ATP8, ND4, and Cyb genes, respectively; four somatic mutations (m.4146A>C, m.4789G>A, m.13345G>A, and m.13498G>A) in ND1, ND2, and ND5 genes, respectively; and one unknown mutation (m.9349T>C) in the COXIII gene, because neither blood nor normal adrenal medulla tissue was available ( Figure 2). These variants are novel and not reported in the MITOMAP database, and in silico analysis predicted them as damaging (Table 2).
Furthermore, among the two novel mutations in tRNA genes, only the m.5658T>C identified in MT-TN (tRNAAsn) was predicted to be pathogenic with the PON-mt-tRNA program (Table 2). Indeed, this m.5658T>C change is located in position 72 of the acceptor stem, which changes the A:U Watson-Crick base pair to a G:U wobble pair. Moreover, the alignment of tRNAAsn sequences from different species showed that the A nucleotide is conserved through different species (Supplementary Figure S3B).

Mitochondrial DNA Depletion and Large Deletions
The mtDNA contents of the 77 PCCs/PGLs and normal adrenal medulla tissues (only available from the French cohort) were measured by ddPCR. A reduction of copy number with more than 50.7% mtDNA content was considered as a depletion. A cutoff value distinguishing between tumors and controls was identified by using ROC curve with 93.8% sensitivity and 89% specificity ( Figure 3A). We found that only 10 tumors (five from the French cohort (FR3, FR30, FR33, F36, and F39) and five from the Swedish cohort (PH40, PH57, PH64, PH65, and PH68)) presented a similar copy number as the normal adrenal medulla tissues. However, the mean mtDNA content of the rest of the tumors was nearly four times lower (3.8-fold) than that in normal PCCs/PGLs tissues (p < 0.0001) ( Figure 3B) (Supplementary Figure S4A).
To assess the mtDNA integrity and mitochondrial rearrangements, we amplified a long fragment (~14 Kbp), including the major arc of mtDNA of the 77 samples and their corresponding controls (Supplementary Figure S4B). The result from the fragment analysis showed that, whilst control samples produced a detectable band of the expected size, 26% (20/77) of PCCs/PGLs presented heteroplasmic somatic deletions with different lengths ( Figure 3C). These deletions are found in the major arc of mtDNA (between the origins of heavy and light strand replication sites), but not in the minor arc (data not shown). The deletions result in the elimination of several genes located in this region, such as tRNA genes and genes encoding respiratory chain proteins. By the amplification of the MT-ND6 gene, we noticed the presence of a heteroplasmic somatic deletion (300 bp) in 39.1% (9/23) of the Swedish tumors, but not in the French cohort or in the controls ( Figure 3D).

Dysregulation of Mitochondrial Genes and Proteins Expression
Gene-expression analysis of the mitochondrial OXPHOS genes showed that there were no differences in mtRNA levels in the tumors presenting novel mitochondrial mutations, as compared to the rest of PCCs/PGLs, in the corresponding mutated genes. On the other hand, analyzed tumors from both the Swedish and French cohorts presented a lower MT-ND2 and MT-ND6 gene expression than normal tissues (p < 0.0001) ( Figure 4A). In addition, the Swedish cohort had a significantly lower gene expression of MT-ND6 when compared to the French cohort (p < 0.05).  Mass spectrometry analysis of isolated mitochondrial proteins showed a slight increase of expression levels of mtOXPHOS proteins (ND4, ND5, COXII, ATP6, and ATP8) in the tumors, presenting different variations in comparison to normal adrenal medulla tissues (p > 0.05) ( Table 3).  Mass spectrometry analysis of isolated mitochondrial proteins showed a slight increase of expression levels of mtOXPHOS proteins (ND4, ND5, COXII, ATP6, and ATP8) in the tumors, presenting different variations in comparison to normal adrenal medulla tissues (p > 0.05) ( Table 3).
A set of protein-expression data across samples were analyzed by using the ProtExAc program. After performing statistical analysis and filtering, we sorted out the differentially expressed proteins and performed enrichment analysis to identify top-scored pathways. A total of 253 proteins were further clustered into 117 groups basis on their enrichment score. The pathway-enrichment analysis showed that most of the proteins were involved in phagosome, oxidative phosphorylation, TCA cycle, pyruvate metabolism, and glutathione metabolism, in addition to other cell-signaling pathways ( Figure 4B).

Absence of 6mA Methylation
Our strategy was to check whether 6mA at the promotor region of bothND2 and ND6 genes could affect the mitochondrial transcription and the low mtDNA copy number. We investigated three sites within previously identified 6mA (Dloop, ND2, and ND6). The 6mA restriction enzyme and semi quantitative PCR identified none of these sites to be methylated (Supplementary Figure S5).

Absence of Mutations in Mitochondrial DNA Integrity Genes
The identification of mtDNA alterations (mtDNA deletions/mtDNA depletion) may also be due to Mendelian disorder, since the maintenance of mtDNA integrity is dependent on several nuclear-encoded proteins. In this context, our microarray data and recent studies (van Gisbergen et al., 2015; Sharma and Sampath, 2019) were used to select the susceptibility nuclear genes, POLG1, C10ORF2, DGUOK, TFAM, and MPV17, which are critically involved in maintaining mtDNA integrity for mutation and gene-expression analysis in our cohorts. No pathogenic mutations were identified, but three common polymorphisms were detected: rs3087374 and rs2307441 in POLG1 in eight and two tumors, respectively; and rs74874677 in DGUOK in two different tumors.

Dysregulation of Mitochondrial Biogenesis and Mitophagy
To assess mitochondrial biogenesis, PCG1α, NRF1, and TFAM gene expression was investigated. Tumors had a very low PCG1α, NRF1, and TFAM gene expression (p = 0.02; p < 0.001 and p < 0.001, respectively) as compared to normal tissues ( Figure 5A). In addition, we investigated two mitophagy markers and we showed a higher expression of LC3a (p = 0.0014) and a lower expression of P62 (p < 0.0001) in those tumors ( Figure 5B).

Discussion
In the present study, we designed a targeted NGS assay for mutation analysis PCCs/PGLs susceptibility genes, followed by the mRNA expression analysis, using a list that distinguishes between pseudohypoxic and RTK/RAS-driven PCCs/PGL mtOXPHOS. In general, we observed that the mtOXPHOS mRNA levels were up lated across the Pheo-Type-based gene-expression signatures, with downregulat MT-ND2/ND6 genes, encoding two subunits of complex I. The differences in com gene expression have been observed previously in many cancers (e.g., lung, colon bladder) [41,42]. Possible explanations are post-transcriptional regulatory mecha and/or differences in mRNA stability among mtDNA-encoded genes. Additionall clear gene expression and transcriptional elements of OXPHOS subunits and ass factors could be factors for the observed changes in the transcript levels [41][42][43]. Th ulation of complex I subunit mRNA is more complex, rather than being a global co lation [43,44], indicating probably an impairment of the OXPHOS complex I. This vation is further substantiated by protein-enrichment analysis, where the OXPHOS way is on the top of the enriched KEGG pathways.
On the other hand, the observed transcriptomic differences between the Swedis French tumors may be explained by the different geographic (ethnic) origins, as wa viously observed between Scandinavian and French PCCs/PGLs during a genetic an of SNPs and PCC/PGL susceptibility [25]. The mitochondrial haplogroup analysi firmed these population differences (Supplementary Figure S2A).
In addition, we performed sequencing of the entire mitochondrial DNA fro PCCs/PGLs patients. We identified eight known pathogenic mutations previously a ated with cancers and 21 novel mutations, where 12 are considered pathogenic. Ov past decades, somatic and germline alterations in nuclear genes and epigenetic ch have been analyzed to identify the molecular basis for tumor formation. Recently, ch

Discussion
In the present study, we designed a targeted NGS assay for mutation analysis of 26 PCCs/PGLs susceptibility genes, followed by the mRNA expression analysis, using a gene list that distinguishes between pseudohypoxic and RTK/RAS-driven PCCs/PGLs and mtOXPHOS. In general, we observed that the mtOXPHOS mRNA levels were upregulated across the Pheo-Type-based gene-expression signatures, with downregulation of MT-ND2/ND6 genes, encoding two subunits of complex I. The differences in complex I gene expression have been observed previously in many cancers (e.g., lung, colon, and bladder) [41,42]. Possible explanations are post-transcriptional regulatory mechanisms and/or differences in mRNA stability among mtDNA-encoded genes. Additionally, nuclear gene expression and transcriptional elements of OXPHOS subunits and assembly factors could be factors for the observed changes in the transcript levels [41][42][43]. The regulation of complex I subunit mRNA is more complex, rather than being a global co-regulation [43,44], indicating probably an impairment of the OXPHOS complex I. This observation is further substantiated by protein-enrichment analysis, where the OXPHOS pathway is on the top of the enriched KEGG pathways.
On the other hand, the observed transcriptomic differences between the Swedish and French tumors may be explained by the different geographic (ethnic) origins, as was previously observed between Scandinavian and French PCCs/PGLs during a genetic analysis of SNPs and PCC/PGL susceptibility [25]. The mitochondrial haplogroup analysis confirmed these population differences (Supplementary Figure S2A).
In addition, we performed sequencing of the entire mitochondrial DNA from 77 PCCs/PGLs patients. We identified eight known pathogenic mutations previously associated with cancers and 21 novel mutations, where 12 are considered pathogenic. Over the past decades, somatic and germline alterations in nuclear genes and epigenetic changes have been analyzed to identify the molecular basis for tumor formation. Recently, changes in mitochondrial cellular content and mtDNA mutations have been proposed as new molecular markers for cancer detection and surveillance [45]. The frequency of mitochondrial mutations in solid tumors is less studied, including bladder and head and neck cancer where hot-spot genes for mtDNA mutations are identified: ND3, ND4, COXI/II/III, 16S rRNA, and deletions in the Dloop region [46,47].
Frameshift mtDNA mutations are rare and were previously reported in six mitochondrial genes: MT-CYB, MT-ND1//4/5/6, and MT-COXIII, with only two frameshift mutations (14510delA and 13384insC) reported in the MT-ND6 gene. Previous studies reported that MT-ND6 frameshift mutations in cell lines manifest a complex I enzymatic deficiency and reduced levels of the assembled complex I, suggesting increased cell death and disease pathogenesis [48]. Using in silico prediction and referring to the functional studies in MT-ND6 gene, we see that the novel 14603insT (Ser24Tyrfs) mutation is considered as pathogenic as the two cited mutations. Besides, screening the MT-NDx genes displayed six somatic and one germline pathogenic mutations, including a nonsense somatic mutation in MT-ND1 gene m.3563G>A (Trp86X). Although the microarray and mass spectrometry data failed to confirm expression alterations, the mutations are predicted to disturb complex I formation and formation of the proton translocation module (P module), as indicated by the Mitomap2 database. Previously, mutations in MT-NDx genes have been associated with different human cancers (e.g., colorectal carcinomas, pancreatic cancer, and oral squamous cell carcinoma) [49]. The MT-NDx mutations inhibit oxidative phosphorylation, increase ROS production, and potentially stimulate the metastatic potential in myeloid leukemia and contribute to carcinogenesis [50].
The screening of MT-COX genes revealed one novel somatic heteroplasmic mutation in MT-COXIII m.9553G>A (Trp116X). This pathogenic mutation results in a truncated protein with the loss of 145aa of the COX3 polypeptide. Mutations in mitochondrial genes encoding MT-COXI/II/III subunits (complex IV) have been associated with numerous cancers. COXIII mutations have mainly been associated with colon and thyroid cancer (https://www.mitomap.org/foswiki/bin/view/MITOMAP/MutationsSomatic) (accessed on 23 September 2021) and suggested to cause an overproduction of superoxide anions by the mitochondrial respiratory chain at COX deficiency [51].
Pathogenicity of the m.5658T>C germline mutation in the tRNAAsn gene has been assessed, suggesting the absence of essential recognition nucleotides. A previous study of the steady-state aminoacylation kinetics of acceptor stem mutant tRNAs demonstrated that replacing any base pair in the acceptor helix with another Watson-Crick base pair affects aminoacylation (or tRNA "charging") efficiency and can elicit mistranslation or lead to a defective tRNA [52].
Certain mutations in mitochondrial DNA produce oxidative phosphorylation defects (preferentially in complex I) that may have pro-tumorigenic effects [53]. For example, MT-ND6 mutations lead to deficient complex I activity and high ROS generation, making the cells highly metastatic. Ishikawa et al., found that replacement of endogenous mtDNA in a poorly metastatic mouse tumor cell line with mtDNA from a highly metastatic cell line (Lewis lung carcinoma) containing two pathogenic mutations in MT-ND6 made the recipient tumor cell acquire the metastatic potential of the transferred mtDNA. Furthermore, these mitochondrial mutations lead to a deficiency in respiratory complex I, associated with overproduction of ROS. Pretreatment of the highly metastatic tumor cells with ROS scavengers suppressed their metastatic potential in mice [54].
The mitochondrial deletions identified in our study could be one of the mitochondrial mutations identified to be causative of human diseases and one among 200 documented pathogenic mtDNA mutations (http://www.mitomap.org/) (accessed on 23 September 2021) [55]. Single and large-scale mtDNA deletions are sporadic and occur because of abnormalities in mtDNA replication and can also arise during the repair of mtDNA damage [56,57]. Large deletions, observed in 16/23 of Swedish tumors and 4/54 of French tumors, result in the loss of several protein-encoding genes and some tRNA genes. Moreover, a deletion of 300 bp was detected only in the Swedish tumors and may have a tumorigenic role similar to the common 4977 bp deletion in the major arc, which has the potential to be a biomarker for cancer occurrence in the tissue [56]. This is interesting, but due to a small sample size, this finding needs to be confirmed in larger PCC/PGL cohorts. Thus, they may affect the mitochondrial polypeptide synthesis and contribute to the impairment of the oxidative phosphorylation and energy metabolism in the respiratory chain [58]. This finding, is in accordance with Neuhaus et al. (2014Neuhaus et al. ( , 2017 observations, where they have showed that catecholamine metabolism could induce mtDNA deletions, causing mitochondrial dysfunction in adrenal medulla and cortex [15,16] Another alteration that may affect cellular function is the copy number and depletion of mtDNA that are observed in PCCs/PGLs. Changes in the tumor mtDNA copy number have been noted in human cancers [59] and found to vary between different cancer forms, being elevated in, for example, primary tumors of head and neck squamous cell carcinoma [60] and in papillary thyroid carcinomas [61], and reduced in breast tumors [62] and gastric cancers [63], relative to normal controls. It has been demonstrated that the reduction of mitochondria leads to tumorigenesis by inducing changes in redox status, membrane potential, ATP levels, gene expression, nucleotide pools, and increased chromosomal instability [64]. Still, it is not yet known whether the depletion of mtDNA is a consequence of losses of mitochondrial function typical of cancer cells, or if it leads to a mutant phenotype adapted to the needs of the cancer cell [65].Further studies are necessary to clarify these issues. On the other hand, our data showed a significant decrease of PCG1α, NRF1, and TFAM gene expression in the analyzed tumors compared to normal tissues. PGC1α is a central regulator of certain nucleus-encoded mitochondrial genes through interactions with multiple transcription factors, including nuclear respiratory factors such as NRF1. NRF1 binds to the promoter and regulate the expression of TFAM and mitochondrial matrix proteins (involved in the replication and transcription of mtDNA) [66,67]. By activating TFAM, PGC-1a could regulate mitochondrial biogenesis, including an increase of mitochondrial copy numbers, activation of the respiratory chain, and augmentation of OXPHOS capacity [68,69]. Our results showed that the reduced number of mitochondria was accompanied by a decrease of the mitochondrial biogenesis, suggested in many recent studies to have a crucial role in the decrease in tumorigenesis or tumor progression [70,71].
Mitochondrial biogenesis, dynamics (fusion and fission), and autophagy (mitophagy) control content and quality of the mitochondria [72]. The impairment of these mechanisms has been associated to cancer, i.e., mitochondrial dynamics [73]. Indeed, DRP1 activation increases the mitochondrial fission in several types of cancers [74] and promotes mitochondrial fragmentation facilitating cancer-cell migration and invasion [75]. Here, we showed a significant overexpression of DRP1 protein (Dynamin-1-like protein) in the analyzed tumors (Table 3), accompanied by decreased mitochondrial biogenesis markers in PCCs/PGLs and a reduced number of mtDNA, which are indicative of mitophagy. To monitor mitophagy, we have extended the expression information on additional biomarkers. Our gene-expression data showed a significant upregulation of LC3a and downregulation of p62. Among the autophagy-related proteins, LC3a is involved in autophagosome formation, and p62 serves as a selective autophagy substrate; both are widely used autophagy markers for monitoring autophagy activity. The complex regulation of LC3a and p62 is suggested to be related to autophagy activity in cancers [76]. These results could suggest the high level of mitophagy [77,78]. We believe that PCCs/PGLs tumor cells may adopt an overactive process of mitochondrial turnover, leading to the selective elimination of dysfunctional mitochondria through mitophagy [79].

Conclusions
Taken together, our results showed that the PCCs/PGLs are caused by germline or somatic mutations in known susceptibility genes in about 50%, and complementary mutations in mitochondrial genes are present in about 17%. These tumors harbored different mitochondrial mutations, deletions, and a reduced copy number and displayed different gene expression patterns, several of which appear without any identified mutation in susceptibility genes, indicating a complementarity and a potentially contributing role in the tumorigenesis PCCs/PGLs.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14020269/s1. Figure S1: (A) Heatmap with top 50 differentially expressed genes in Swedish and French PCCs/PGLs after gene-set enrichment analysis (GSEA). (B) GSEA enrichment plots of canonical pathways in Swedish tumors: Reactome metabolism of proteins and Reactome GPCR downstream signaling. Figure S2: Haplogroup distribution in all studied PCCs/PGLs tumors (our study) (A) compared with a healthy control cohort (Reference group) and other types of cancer (B). Figure S3: Figure S5: Semi-qPCR results of 4 selected CATG sites. After DpnI digestion, semi-qPCR was performed by using specific primers covering these sites. The ratio was calculated between the PCR product concentration from digested and undigested DNA samples (mean ± standard error). Table S1: Clinical information of patients with pheochromocytomas and paragangliomas from Nancy, France, and from Linköping, Sweden. Table S2: List of genes selected for the gene panel sequencing. Table S3: cDNA and DNA primers used for nuclear genes amplification, Sanger sequencing, and RT-qPCR. Table S4: Statistics of the GSEA of Swedish tumors compared to French tumors. Enrichment pattern with FDR lower than 25% was considered significant. Table S5: Association between mtDNA haplogroups bins and PCCs/PGLs.  Informed Consent Statement: Informed consent was obtained from all subjects involved in this study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. buffer containing 20 mM Hepes-NaOH, pH 7.5, 10 mM KCl, 1.5 mM MgCl 2 , 1 mM DTT, 250 mM sucrose, and protease inhibitor (Halt Protease Inhibitor Cocktail, Thermo-Fisher); homogenized (TissueLyser (Qiagen) operated for 2 min at 20-30 Hz, two times); and centrifuged at 1000× g for 15 min, at 4 • C. The pellet was discarded, and the supernatant was centrifuged at 12,000× g for 15 min, at 4 • C. The resulting pellet containing mitochondria was resuspended in 50 mM Tris buffer, pH 7.4, with protease inhibitor, and the protein concentration was determined by using a BCA Protein assay kit (Thermo-Fiher).

Appendix B. Haplogroups Analysis
Whole-mitochondrial sequencing allows us to perform the mitochondrial haplogroup analysis via MitoTool database (http://www.mitotool.org/genome.html) (accessed on 23 September 2021). Haplogroups are mtDNA-specific sequence polymorphism variations. We compared our result by using the classification published in the previous report [80]. In addition, due to the rarity of some haplogroups in our cohort and the small sample sizes in some subgroups, we combined mitochondrial haplogroups into larger bins based on the haplogroup evolutionary tree, including HV (H+V+HV), JT (J+T), and UK (U+K) for comparison of our cases and controls from the literature [81].

Appendix C. Sanger Sequencing
To confirm all mutations and their status in corresponding normal tissue, PCR products were amplified and sequenced on both strands, using specific primers (Supplementary  Table S1) [22,82]. Direct sequencing of PCR products was performed with the ABI Prism BigDye Terminator 3.1 Cycle Sequencing Ready Reaction Kit (ABI PRISM/Biosystems), and the products were resolved on ABI 3500 sequencer (Applied Biosystems). The blast homology searches were performed by using the program available at the National Center for Biotechnology Information website in comparison with the updated consensus Cambridge sequence.

Appendix D. In Silico Mutational Analysis
The amino acid sequence alignment for evolutionarily conserved status was analyzed by using the ClustalW program (http://www.ebi.ac.uk/Tools/msa/clustalw2) (accessed on 24 September 2021), and corresponding sequences of tRNA gene in several species were obtained from the Mammalian Mitochondrial tRNA Genes database (http://mamit-trna.ustrasbg.fr/) (accessed on 24 September 2021).
Performances of in silico prediction tools for non-synonymous mtDNA variants were made according to Mitomap and assessed with 19 different prediction tools gathered in MitImpact2 (http://mitimpact.css-mendel.it/) (accessed on 24 September 2021), including the known tools PolyPhen, Sift, PROVEAN, PhD-SNP, and Mutation Taster, as well as recently developed tools for mtDNA as MToolBox [83], the meta-predictor APOGEE [84], or Mitoclass.1 [85]. Few tools are dedicated to mitochondrial tRNAs; for example, PON-mt-RNA (http://structure.bmc.lu.se/PON-mt-tRNA/) (accessed on 24 September 2021) is a multifactorial score associating 12 features, including evolutionary conservation; primaryto-tertiary structures; and functional assays, including biochemistry and histochemistry [86]. The interpretation of ribosomal variants is based on combined conservation information with structural data, using mFold (http://unafold.rna.albany.edu/?q=mfold/RNA-Folding-Form) (accessed on 24 September 2021 [87]. To predict the effect of the mutation on the 2D structure, we used the PROTTER software V1.0, an open-source tool for visualization of proteoforms and interactive integration of annotated and predicted sequence features, together with experimental proteomic evidence (http://wlab.ethz.ch/protter/) (accessed on 24 September 2021). Using cBioportal and based on the genomic alterations (mutations and mitochondrial copy-number alterations), the portal automatically generates a report summarizing genomic data across all patients through an OncoPrint [88,89].