Genome-wide association study of mitochondrial copy number

Abstract

Mitochondrial DNA copy number (mtDNAcn) variation has been associated with increased risk of several human diseases in epidemiological studies. The quantification of mtDNAcn performed with real-time PCR is currently considered the de facto standard among several techniques. However, the heterogeneity of the laboratory methods (DNA extraction, storage, processing) used could give rise to results that are difficult to compare and reproduce across different studies. Several lines of evidence suggest that mtDNAcn is influenced by nuclear and mitochondrial genetic variability, however this relation is largely unexplored. The aim of this work was to elucidate the genetic basis of mtDNAcn variation. We performed a genome-wide association study (GWAS) of mtDNAcn in 6836 subjects from the ESTHER prospective cohort, and included, as replication set, the summary statistics of a GWAS that used 295 150 participants from the UK Biobank. We observed two novel associations with mtDNAcn variation on chromosome 19 (rs117176661), and 12 (rs7136238) that reached statistical significance at the genome-wide level. A polygenic score that we called mitoscore including all known single nucleotide polymorphisms explained 1.11% of the variation of mtDNAcn (p = 5.93 × 10⁻⁷). In conclusion, we performed a GWAS on mtDNAcn, adding to the evidence of the genetic background of this trait.

Introduction

Mitochondria are organelles involved in the regulation of critical cellular functions such as energy production via the oxidative phosphorylation reaction, apoptosis and calcium homeostasis and are responsible for the production of reactive oxygen species. Mitochondria possess their own genome (mtDNA), which is maternally inherited. mtDNA content has been shown to be a proxy for mitochondrial status and function, and is consequently an attractive biomarker due to its relative ease of measurement (1).

In each eukaryotic cell there can be hundreds or thousands of copies of their genomes. Mitochondrial DNA copy number (mtDNAcn) can be readily measured by quantitative real-time PCR (qPCR) using DNA extracted from peripheral blood or other tissues, although in practice most measurements are performed on material from peripheral blood due to ease of accessibility. mtDNAcn varies by cell type, but there is generally a correlation in the amount of mtDNA in different cell types. Therefore, mtDNAcn measured in circulating leukocytes could represent a good indicator of the average amount of mtDNAcn in other tissues (2).

Many studies investigated the possible role of mtDNAcn as a risk factor for several cancer types (3), including breast (2,4), colorectal (5), lung (6), prostate (7), pancreatic cancers (8,9) and B-cell lymphomas (10). Besides cancer, mtDNAcn has been associated with a wide range of phenotypes, such as myocardial function (11), cardiovascular disease (12), male infertility (13), depression (14), posttraumatic stress disorder (15), chronic stress and stress-related health behaviours (16). In general, since mitochondria are at the crossroads of many fundamental cellular processes, it can be expected that mtDNAcn might be associated with a wide range of physiological and pathological phenotypes.

mtDNAcn is influenced by exogenous as well as endogenous factors, including genetic factors since it appears to have a high heritability, estimated at 65% (17,18). Single nucleotide polymorphisms (SNPs) in candidate gene and genome-wide association studies (GWASs) were reported to be associated with mtDNAcn (17,19–24). Curran and colleagues in 2007 performed the first GWAS on mitochondrial content. In their analyses 1259 Mexican American individuals were included, and mtDNAcn was estimated using qPCR (17). López and colleagues performed the second GWAS on mitochondrial content measured through qPCR using 21 Spanish families consisting of 386 individuals spanning across five generations. They observed several suggestive loci that did not reach genome-wide significance (21). Cai and colleagues used low-coverage sequence data from 5224 women with major depressive disorder and 5218 healthy controls to compute the normalised numbers of reads mapping to the mitochondrial genome as a proxy for the amount of mtDNA. They identified two SNPs (rs11006126 and rs445) associated with mtDNAcn (22). Workalemahu and colleagues evaluated, in 471 Peruvian women, the association between mtDNAcn measured by qPCR and 119 629 SNPs (25). A small scale GWAS on mtDNA performed on 301 subjects from the Chinese population using qPCR identified a few suggestive loci, but none reached genome-wide significance (p < 5x10⁻⁸) (20). Yamamoto and colleagues analysed mtDNA using the whole genome sequence of 1928 Japanese subjects and did not report any associations (26). In addition, Guyatt and colleagues performed a GWAS on very heterogeneous populations of European descent using different setups of qPCR, suggesting novel candidates which did not reach the genome-wide significance (23). Finally, Hägg and colleagues conducted a study using data from UK Biobank obtaining the mtDNAcn from genotyping microarray probe intensities and reported 50 novel genome-wide significant loci (24).

The majority of these studies analysed together more than one population group, or small groups of subjects collected in different context or included only females, making thus the studies less homogeneous among them. These factors make it difficult to compare different studies and provide a global interpretation of the results. In addition, the application of different mtDNAcn quantification methods contributes to increasing the heterogeneity between studies (27,28). qPCR is the current de facto standard and has been the most widely used method for measuring mtDNAcn. More recently, new technologies such as droplet digital PCR technology, the use of genotyping assay intensity or the numbers of reads obtained from whole exome sequencing and whole-genome sequencing data have become available (28,29).

In addition to the different study designs and quantification techniques, other factors could increase the heterogeneity of the results between different studies, such as the nature of the samples from which the DNA is extracted, the extraction techniques used, or the content of the different cells in the samples (white and red blood cell percentage, and platelet count). (11,28). The absence of uniformity in sample handling, quantification techniques and study design has produced contrasting results (30,31). Indeed, both positive and negative associations have been found for mtDNAcn and aging (32,33), type 2 diabetes (34,35), pancreatic cancer (8,9), and breast cancer (4,6,36). For example, retrospective studies are prone to reverse causation bias if the disease or its treatment influence the mitochondrial count.

To increase our knowledge on the genetics of mtDNAcn, we conducted a study on a homogeneous population using the same measurement technique. We performed a GWAS on 6836 unrelated German individuals from the ‘Epidemiologische Studie zu Chancen der Verhütung, Früherkennung und optimierten THerapie chronischer ERkrankungen in der älteren Bevölkerung’ (ESTHER) prospective cohort study (37), to identify novel genetic variants associated with mtDNAcn and to compute a polygenic score (PGS) which we called ‘mitoscore’.

Results

Associations between SNPs and mtDNA copy number

In the ESTHER population we observed two genome-wide statistically significant associations with mtDNAcn: one on chromosome 13q12.2 (rs58413920) and one on chromosome 9q31.1 (rs10124625). Both minor alleles, A for rs58413920 and T for rs10124625, were associated with increased mtDNAcn: beta = 0.259; 95%CI: 0.176–0.343; p = 1.30 × 10⁻⁹ and beta = 0.072; 95% CI: 0.046–0.097; p = 3.36 × 10⁻⁸ respectively. rs58413920 is situated near the small nuclear RNA (snRNA) pseudogene RNU6-70P (8496 bp downstream of the SNP), while rs10124625 lies between the 5S ribosomal pseudogene RNA5SP291 (122 831 bp upstream the SNP) and the uncharacterised lncRNA LOC101928523 (187 347 bp downstream of the SNP). All the results, alongside information on the nearest genes of each SNP, chromosome position and allele frequency, are shown in Supplementary Material, Table S1. A Manhattan plot is shown in Figure 1, and Supplementary Material, Figure S1 shows regional plots for the two new genome-wide significant loci. In this study were included as replication set the summary statistic results of the recent GWAS by Hägg and colleagues performed using 295 150 participants from the UK Biobank (24). We performed a meta-analysis of the results obtained in the ESTHER and UK-Biobank studies excluding the SNPs already known to be associated with mtDNAcn (and those associated with them r² > 0.8) and keeping only the SNPs with p < 0.05 in the ESTHER study. The 13q12.2 (rs58413920) SNP (MAF = 0.01 in Europeans) and its proxies were not present in the summary statistics of the UK-Biobank GWAS (24) due to the MAF threshold (0.05) set by Hägg and colleagues, instead, 9q31.1 (rs10124625) was present and with beta = 0.007 and p = 0.045.

Figure 1

Manhattan plot of the whole genome association analysis of mtDNAcn.

In the meta-analysis we identified two novel SNPs that reached genome-wide p-value threshold: 12q15-rs7136238 (beta = −0.013; p = 3.83 × 10⁻⁸), and 19p13.11-rs117176661 (beta = 0.034; p = 2.73 × 10⁻⁹). The association results are reported in Table 1.

Table 1

Associations between SNPs and mtDNAcn

							ESTHER				UK-Biobank				Meta-analysis
CHR	BPa	SNP	EAb	NEA	MAFc  (ESTHER)	MAFc  (EUR)	BETAd	SEd	Pd	N	BETA	SE	P	N	BETA	SE	P
12	69 752 509	rs7136238	T	G	44%	48%	−0.020	0.010	0.041	6735	−0.013	0.003	2.68 × 10−7	286 381	−0.014	0.003	3.83 × 10−8
19	17 445 208	rs117176661	T	C	4%	6%	0.090	0.026	4.71 × 10−4	6623	0.031	0.006	1.13 × 10−7	286 381	0.034	0.006	2.73 × 10−9

^avariant position based on Genome Reference Consortium Human Build 37 (GRCh37); ^bEA: effect allele; ^cMAF: minor allele frequency; ^dLinear regression models to examine the association between SNPs and mtDNAcn were adjusted by sex, BMI, fruit intake frequency, pack-years of cigarette smoking.

The two polymorphisms identified in the meta-analysis were independent (r² < 0.15 in EUR subjects of the 1000 Genomes project) with SNPs reported in other studies. This independence was confirmed from the results obtained from the functional mapping and annotation of genome-wide association studies (FUMA) portal. The results obtained from FUMA reported three independent genome-wide lead SNPs in the 1 000 000 bp region around rs117176661 (rs35006574; rs67397200; rs117176661) on chromosome 19, and one around rs7136238 (rs7136238) on chromosome 12. FUMA results on lead SNPs are reported in Supplementary Material, Table S2. Both variants identified through the meta-analysis are reported in GTEx as expression quantitative trait loci (eQTLs) for six genes (three for each SNP), in a large number of human tissues. The detailed results from bioinformatic annotation tools are reported in Table 2. RegulomeDB assigns a score of 4 and 5 to 19p13.11-rs117176661 and 12q15-rs7136238, respectively, indicating a moderate chance for the SNPs to have a functional role.

Table 2

Results of bioinformatics analyses

CHR	BPa	SNP	Geneb	Rdb-Rank	eQTLd	eQTL-tissue	Gene Name	Distance bp	DnIfe
12	69 752 509	rs7136238	LYZ [UG] 4495	5	LYZ	33 tissues
			YEATS4 [DG] 974		RP11-1143G9.5	14 tissues
					YEATS4	48 tissues
							CCT2; MIR3913–1; MIR3913–2	225 000	3.02
19	17 445 208	rs117176661	ANO8 [OG]	4	ANKLE1	6 tissues
					MRPL34	4 tissues	MRPL34	30 000	3.68
					OCEL1	4 tissues
							B3GNT3	460 000	8.62
							INSL3	485 000	5.47
							NR2F6	90 000	4.1
							DDA1	25 000	3.97
							FAM129C	185 000	3.78
							MAP1S	385 000	3.77
							BST2	70 000	3.76
							GLT25D1	220 000	3.67
							F2RL3	450 000	3.55
							CPAMD8	310 000	3.51
							NXNL1	125 000	3.08

^aPosition: variant position based on Genome Reference Consortium Human Build 37 (GRCh37). ^b [OG] overlapped gene; [UG] nearest upstream gene: [DG] nearest downstream gene; distance from nearest gene; information obtained from SNPnexus. ^c RegulomeDB rank. ^d GTEx database Release V8 (data obtained from the GTEx Portal on 01 March 2021). ^e DIV3 -Distance-normalised interaction frequency.

On chromosome 19, we reported a novel variant associated with mtDNAcn. Near this locus, Hägg and colleagues reported a variant that reached the genome-wide significance level (rs10419397, 51 955 bp distant from rs117176661), while Guyatt and colleagues identified another variant (rs10424198, 33 612 bp distant from rs117176661) lying between these two SNPs. Although independent of each other (rs10419397 r² = 0.014 D′ = 0.92; rs10424198 r² = 0.017 D′ = 1.00 with rs117176661 in EUR), these SNPs are physically very close, indicating the regulatory potential of this genomic region on chromosome 19. We used the 3D Interaction Viewer and Database (3DIV) (38) to obtain epigenetic and spatial chromatin interaction annotation of the region between those three polymorphisms. For all the analyses, we considered a window of 500 000 bp and a distance normalized interaction frequency (DNIF) threshold of 3 using the data obtained from IMR90 fibroblast cell line. 3DIV indicates that rs117176661 belongs to a different regulatory topologically associating domain compared to rs10419397 and rs10424198. 19p13.11-rs117176661 interacts with 12 genes. In contrast, the variant identified in the UK-Biobank study (rs10419397) interacts with just one gene when analysed with the same parameters. 12q15-rs7136238 interacts with chaperonin containing TCP1 subunit 2 (CCT2) gene and two micro-RNA (MIR3913–1, MIR3913–2). The interaction values and the list of genes are reported in Table 2.

We computed in the ESTHER GWAS a PGS that we called mitoscore using 44 SNPs (out of the 50 loci reported by Hägg as genome-wide significant, six were not present in the imputed dataset). The mitoscore explains 1.05% (Nagelkerke pseudo-r² = 0.01049; p = 3.56 × 10⁻⁶) of the variation of mtDNAcn, while the PGS generated using the 44 variants plus the two novel variants identified in this study explains 1.11% (Nagelkerke pseudo-r² = 0.01107; p = 5.93 × 10⁻⁷) of the variation. The mitoscore (44 SNP) was categorized in quintiles, and with a logistic regression the fifth quintile (high mtDNAcn) was compared with the first quintile (low mtDNAcn) showing a statistically significant association (beta = 6.03; 95% CI = 3.87–8.19; p = 4.33 × 10⁻⁸) with the mtDNAcn. The results of the analyses of mitoscore are reported in Supplementary Material, Table S3.

Discussion

The variation in mtDNAcn measured in circulating blood has been associated with various human traits and diseases (2,4–11,13–16). In the vast majority of the reports, mtDNAcn has been measured through qPCR quantification, which is the de facto standard for relative quantification of mtDNAcn (28). There is large heterogeneity in the reported results, which could be largely due to the different study designs used and the differences in DNA extraction protocols, storage and processing (30,31,39). We performed an additional GWAS considering mtDNAcn as the outcome and to validate our results, we used the summary statistics results of the study performed on 295 150 participants from the UK Biobank. We have identified two novel independent genome-wide significant associations for mtDNAcn.

The first of these associations is between the minor allele (T) of rs117176661 on chromosome 19 and increased mtDNAcn. This SNP is relatively close to two other SNPs already known to be associated with mtDNAcn. However, the bioinformatic analyses indicate that this SNP is not in LD with the known variants and is located in a different regulatory topologically associating domain. Graphical representation of the locus on chromosome 19 and the chromatin interaction is shown in Supplementary Material, Figure S2. The presence of two different regulatory domains tagged by two SNP identified in independent studies suggests the locus’s potential role and at the same time that there is more than one regulatory region inside the locus that can modulate separately and in synergy mtDNAcn phenotype. 19p13.11-rs117176661 has a RegulomeDB rank of 4, indicating that it lies in a DNAse accessible region and probably alters possible transcription factor’s binding sites. It is an eQTL of three genes in several tissues, and one of these is the mitochondrial ribosomal protein L34 (MRPL34) gene. The potential regulatory role of this SNP on MRPL34 is supported by 3DIV tool, which reports a chromatin interaction with the MRPL34 gene region with a DNIF of 3.8. The data obtained from the bioinformatics analyses are potentially in agreement with our results: we observed that 19p13.11-rs117176661-T is associated with an increase in the relative mtDNAcn, and the same allele is associated with increased expression of the MRPL34 gene, which is involved in the synthesis of mitochondrial ribosomal proteins. The increment of the synthesis of an element essential for the mitochondrial function could be potentially correlated with the number of mitochondria, suggesting a potential new nuclear locus directly involved in the regulation of mitochondrial homeostasis.

12q15-rs7136238 has a RegulomeDB rank of 5 indicating a low regulatory role and is a multi-tissue eQTL for two genes, lysozyme (LYZ), YEATS domain containing 4 (YEATS4) and a long noncoding RNA (RP11-1143G9.5). 3DIV tool reports that this SNP interact with CCT2 gene and two micro-RNAs (MIR3913–1, MIR3913–2) located at a distance of 225 000 bp. Although the annotation tools report some evidence for a functional role of 12q15-rs7136238, the available information does not allow us to hypothesize a possible role in regulating the number of mitochondria.

For leukocyte telomere length, a marker for several human diseases that suffers from the same technical problems as mtDNAcn, a genetic score with telomere length has been identified, validated and applied in several epidemiologic studies (40–45). We sought to do the same for mtDNAcn carrying out a homogeneous study. The mitoscore generated with 44 of the 50 loci identified in UK-Biobank study explains 1.05% (p = 3.56 × 10⁻⁶) of the phenotypic variability and 1.11% (p = 5.93 × 10⁻⁷) adding the two new SNP to the PGS. Many studies have reported a small r² value < 5% and in some cases < 1% for various PGS in different traits (46–50), therefore reporting an r² of 1.11% with the first PGS generated for a complex trait as mtDNAcn can be considered a starting point for future studies to implement new variants to the score. Adding the two new SNPs increased the variance (+0.6%) explained by the score, suggesting that the SNPs we identified are relevant predictors of this phenotype.

A potential limitation of this study is the absence of blood cell counts, a factor associated with mtDNAcn as reported by a recent study (24). One of the strengths of this study is that it has been carried out on a large number of subjects of the same ethnicity and of both sexes, without any kind of internal stratification and including several environmental/lifestyle potential confounders (i.e. dietary behaviours). Our study (in the discovery phase) is extremely homogeneous because we used a measurement technique that is considered the de facto standard for these studies, on all the 6836 samples that were collected, stored and extracted in the same way and analysed at the same time.

In addition, another point of strength is the two-step approach (discovery and a replication phase), the latter performed in one of the largest populations ever tested for this phenotype.

In conclusion, we identified two new genetic variants (19p13.11-rs117176661 and 12q15-rs7136238) involved in the variation of mtDNAcn, and we propose for the first time the mitoscore as a tool to estimate the genetic basis of mtDNAcn variation. We envision, in the future, the use of this score as a proxy of direct measurement of mitochondria, reducing time, cost, and possible bias inherent to the direct measurement of mtDNAcn, thus improving reproducibility of mtDNAcn analyses in epidemiological studies.

Table 3

Characteristics of the study population

	Total	Male	Female
Number of Individuals	7640	3462	4178
Median age, years (25th–75th percentile)	62 (57–67)	63 (57–67)	62 (57–67)
50–64 years at baseline	4700 (62%)	2126 (61%)	2574 (62%)
65–75 years at baseline	2940 (38%)	1336 (39%)	1604 (38%)
Number of individuals with information on fruit intake frequency	7381	3352	4029
more than once/day	2198	687	1511
once/day	2415	1095	1320
more than one/week	2205	1220	985
once/week	290	192	98
less than once/week	221	126	95
no	52	32	20
Number of individuals with information on vegetable intake frequency	7323	3337	3986
more than once/day	266	79	187
once/day	1626	531	1095
more than one/week	4614	2263	2351
once/week	681	377	304
less than once/week	120	73	47
no	16	14	2
Number of individuals with information on BMI	7631	3457	4174
Underweight (BMI < 18,5)	34	10	24
Normal range (BMI < 25)	2035	757	1278
Overweight (BMI < 30)	3612	1852	1760
Obese (BMI > = 30)	1950	838	1112
Median BMI (25th–75th percentile)	27.71 (25–30)	27.86 (25–30)	27.58 (24–30)
Number of individuals with information on alcohol consumption	7640	3462	4178
Abstainer	2938	854	2084
Women 0–19.99 g/d or Man 0–39.99 g/d	4211	2359	1852
Women 20–39.99 g/d or Man 40–59.99 g/d	386	182	204
Women > =40 g/d or Man > =60 g/d	105	67	38
Average g of alcohol consumed for a day (25th–75th percentile)	9.00 (0–13)	13.88 (1.70–19.71)	4.95 (0–6.29)
Number of individuals with information on physical activity	7640	3462	4178
Inactive: < 1 h of physical activity/week	1628	537	1091
Low: 1-2 hrs physical activity/week	3496	1531	1965
Medium or high: > = 2 h of vigorous and > = 2 h of light physical activity/week	2516	1394	1122
Number of individuals with information on smoking behavior	7424	3401	4023
non-smoker	3702	1014	2688
former smoker	2465	1719	746
current smoker	1257	668	589
Pack-years of cigarette smoking (25th–75th percentile)a	22.41 (7–34)	25.07 (9–38)	17.65 (4.8–27.5)

^avalues measured on former and current smokers’ groups

Materials and Methods

Study population

The study was based on the (ESTHER) cohort study, whose details have been reported elsewhere (37). Briefly, 9940 individuals, aged 50–75 years, have been enrolled between 2000 and 2002, and followed for 17 years. After providing informed consent, information about lifestyle, smoking habits, diet, and sociodemographic characteristics were collected at baseline as well as at follow-ups. The ESTHER study was approved by the ethics committees of University of Heidelberg and the state medical board of Saarland, Germany.

In this study we included the information on age (reported as continuous variable), sex (dichotomous variable), BMI (continuous variable), smoking status, pack-years of cigarette smoking among current and former smokers (continuous variable), physical activity (PA; classified as 0: ‘inactive’ (reference) < 1 hr PA/week; 1: ‘low’ 1-2 hrs PA/week; 2: ‘moderate’/‘high’ ≥2 hrs PA/week), fruit intake frequency (1: more than once/day; 2: once/day; > 3: more than one/week; 4: once/week; 5: less than once/week; 6: no), vegetable intake frequency (1: more than once/day; 2: once/day; > 3: more than one/week; 4: once/week; 5: less than once/week; 6: no), and alcohol intake (continuous variable, number of grams ingested daily). Main characteristics of the study population are reported in Table 3.

DNA extraction and sample preparation

Blood samples were taken during a routine health examination and stored at −80°C until analysis. DNA from whole blood samples was collected using a salting out procedure. DNA concentration of each sample was measured using PicoGreen dsDNA quantitation assay (Invitrogen|ThermoFisher, Waltham, MA, USA) to intercalate the dsDNA and the Tecan GENios Multireader microplate reader to read the UV absorbance. The measured DNA of each individual was diluted with water to the final concentration of 5 ng/μL.

qPCR measurement of mtDNA copy number

mtDNAcn measurement was performed in 7640 study subjects. Real-time multiplex polymerase chain reaction was used to quantify the copy number of the mitochondrial gene NADH dehydrogenase, subunit 1 (ND1), and of the nuclear single copy gene albumin (ALB) as housekeeping reference. All the DNA samples were diluted to a concentration of 2.5 ng/μL in a 96 well plates. The triplicates of each sample were obtained dispensing with an electronic multichannel pipette 2 μL of the sample in 3 wells of a 384-well plate (5 ng DNA of each individual in each well). We previously reported the details of the quantification technique (8). The PCR was performed using a Viia-7 sequence detection system (Applied Biosystems) to acquire the cycle threshold (Ct) values for copy number of ND1 and ALB gene.

Individual Ct values that deviated from the average of the triplicates by more than 5% of the standard deviation were discarded. This affected 467 (6%) samples. Furthermore, samples whose average Ct was not within the standard curve range were excluded and not included in the analyses (N = 337, 4%). Standard curves were generated in each plate with a serial dilution (1:2) from 30 ng to 0.47 ng of genomic DNA pooled from 50 random individuals belonging to the study. The standard pool (STP) was renewed every 4 ~ 5 plates (using however the same subjects), to preserve the quality of the pooled DNA during the laboratory measurement. For each reaction the RT-PCR efficiency (E) was calculated using standard curve points in the exponential phase according to the equation: E = 10^[−1/slope].

After exclusion due to quality control (QC) failure, the mtDNAcn was obtained for 6836 individuals as the ratio between ND1/ALB copy numbers, based on the calculation introduced by Pfaffl (51). This approach is suited to qPCR results lacking identical efficiency between the amplification reaction of ND1 and ALB, using the Ct of the standard curve to the equivalent of 5 ng of DNA as a calibrator.

Genotyping and imputation

Genome-wide SNP genotyping was performed on all study subjects using the Illumina Infinium OncoArray and Global Screening Array BeadChips (Illumina, San Diego, CA, USA) as previously described in detail (52). Genotyping quality controls followed the protocols reported by Anderson et al. (53). Genotyping data were phased using SHAPEIT2 and imputed using the Michigan Imputation Server, and MiniMac 4 was used to impute to the haplotype reference consortium (HRC) Version r1.1 2016 reference panel (54,55).

Correlation analyses

We performed a correlation analysis (Pearson correlation) to estimate the linear relation between mtDNAcn and the covariates included in our study (age, sex, BMI, physical activity, alcohol consumption, pack-years of cigarettes smoking, fruit, and vegetable intake frequency). The plates used for each qPCR were classified in group on the base of the STP included in them, and this classification was used and evaluated as a covariate. The covariates that showed a statistically significant correlation (p < 0.05) were included in the association analyses as adjustment variables (Table 4).

Table 4

Correlations between mtDNA copy number and selected covariates

Variable	ra	p-valueb
Age (years)	−0.019	1.20 × 10−1
Alcohol consumption (g/day)	−0.012	3.64 × 10−1
BMI	−0.036	2.77 × 10−3
Cigarette smoking (pack-years)	−0.059	3.37 × 10−6
Fruit intake frequency	−0.043	4.81 × 10−4
Physical activity	0.001	9.28 × 10−1
Sex	−0.065	6.75 × 10−8
Vegetable intake frequency	−0.014	2.43 × 10−1
STP (standard pool)	0.013	2.74 × 10−1

^ar: Pearson correlation coefficient; ^bp-values in bold are statistically significant (p < 0.05).

Genome-wide association analysis

PLINK1.9 software was used to assess genotyping QC on the subjects with mtDNAcn available. Population structure was assessed by principal component analysis. SNPs with minor allele frequency (MAF) < 0.01, or evidence for violations of Hardy–Weinberg equilibrium (p < 10⁻⁵) were excluded, leaving for the analysis 6 386 351 SNPs of 6836 individuals. The lambda value, used to estimate the inflation rate, was 0.995 and a QQ plot is reported as Figure 2. We used linear regression models to examine the association between SNPs and mtDNAcn as a quantitative trait, adjusting for the variables associated with the mtDNAcn in the correlation analyses (sex, BMI, fruit intake frequency, pack-years of cigarette smoking), the top eight principal components and the STP. For statistical significance, we used a threshold of p < 5 × 10⁻⁸.

Figure 2

Quantile-quantile (QQ) plot for the whole genome association analysis of mtDNAcn.

Summary statistics of UK biobank GWAS

The summary statistics of the GWAS performed using 295 150 participants from the UK Biobank (24) were used as replication. The mtDNAcn was estimated from the intensities of array genotyping probes on the mitochondrial chromosome. The betas, confidence intervals and p-values were obtained by logistic regression adjusted for the first ten principal components, age at baseline, sex, genotyping batch, genotyping missingness/call rate, and neutrophil percentage lymphocyte percentage and white blood cell count. The details of the study are reported in the original publication (24).

GWAS meta-analysis

A meta-analysis of the results obtained in the ESTHER GWAS and in the UK-Biobank GWAS was performed. The association results (beta, and standard error) for each SNP were combined with the METAL software (v2011-03-25) adopting a fixed effects model. The meta-analysis included 2 775 252 SNPs with genotype data available from the two studies. A p-value of 5 × 10⁻⁸ was used to establish a threshold for genome-wide significance.

Polygenic scores

We computed a PGS to estimate the relative mitochondrial content. The score was generated using the effect sizes (betas) reported in the summary statistics of the GWAS performed on the UK-Biobank population and was tested in the ESTHER population.

The score was calculated using PLINK software, and the association with mtDNAcn was tested with a generalised linear model using the quantification values measured in the ESTHER population.

We generated two scores, one including the variants identified by Hägg and colleagues. Of the 50 genome-wide significant loci reported by Hägg, 44 were present in the ESTHER GWAS and were included in the score. A second score was computed, adding to the 44 variants the new variants identified in this study. The PGS for each study subject was generated by adding up the product of the effect (beta) of each variant by the number of effect alleles. To take into account that not all subjects had a 100% call rate, the PGS was divided by the number of genotypes available for each subject (average subject call rate 0.99). The analyses for all PGSs were adjusted by the same covariates used for the genome-wide analyses of individual SNPs.

Bioinformatic analysis

Bioinformatic tools were used to evaluate the functional relevance of each SNP with statistically significant results in the genome-wide association analysis. Independent signals within each locus were defined with a correlation (r²) of less than 0.1 to the variant with the lowest p-value using the FUMA platform (56). FUMA was also used to integrate the functional annotation of associated variants. SNPnexus was used as database for functional annotation of SNPs (57). Regulome DB was used to identify the regulatory potential of each SNP (58), and Genotype-Tissue Expression GTEx Analysis Release V8 (data obtained from the GTEx Portal on 01 March 2021) was used to identify the relationship between the SNPs and the level of expression in the nearby genes (59). Finally, we used the web-tool 3D-genome interaction viewer (60,61) to evaluate the chromatin interaction and the epigenetic annotation in the newly identified risk loci.

Acknowledgements

We thank Drs. Sara Hägg and Felix Grassmann (Department of Medical Epidemiology and Biostatistics, Karolinska Institutet, Stockholm, Sweden) for kindly sharing with us the summary statistics of the study reported in Hägg et al, 2021.

This work was supported by intramural funding of German Cancer Research Center (DKFZ) and of the University of Pisa. The ESTHER study was supported by grants from the Baden-Württemberg Ministry of Science, Research and Arts, the German Federal Ministry of Education and Research, the German Federal Ministry of Family, Senior Citizens, Women and Youth, the Saarland Ministry of Social Affairs, Health, Women and Family, and the Network Aging Research at Heidelberg University.

Conflict of Interest statement. The authors declare no competing interests.

Author notes

Manuel Gentiluomo and Matteo Giaccherini authors share the first position

Federico Canzian and Daniele Campa authors share the last position