Skip to main content
Download PDF

The mitoepigenome responds to stress, suggesting novel mito-nuclear interactions in vertebrates

BMC Genomics volume 24, Article number: 561 (2023) Cite this article

Abstract

The mitochondria are central in the cellular response to changing environmental conditions resulting from disease states, environmental exposures or normal physiological processes. Although the influences of environmental stressors upon the nuclear epigenome are well characterized, the existence and role of the mitochondrial epigenome remains contentious. Here, by quantifying the mitochondrial epigenomic response of pineal gland cells to circadian stress, we confirm the presence of extensive cytosine methylation within the mitochondrial genome. Furthermore, we identify distinct epigenetically plastic regions (mtDMRs) which vary in cytosinic methylation, primarily in a non CpG context, in response to stress and in a sex-specific manner. Motifs enriched in mtDMRs contain recognition sites for nuclear-derived DNA-binding factors (ATF4, HNF4A) important in the cellular metabolic stress response, which we found to be conserved across diverse vertebrate taxa. Together, these findings suggest a new layer of mito-nuclear interaction in which the nuclear metabolic stress response could alter mitochondrial transcriptional dynamics through the binding of nuclear-derived transcription factors in a methylation-dependent context.

Peer Review reports

Background

Mitochondria are key mediators of the physiological and pathophysiological cellular responses to environmental stressors. The association between mitochondrial dysfunction and neurological disorders as well as cancer, Down syndrome, diabetes and infertility is becoming increasingly recognized [1, 2]. Importantly, mitochondrial dysfunctions are seldom linked to mitochondrial DNA (mtDNA) sequence alterations [1, 3], fuelling a burgeoning interest in the mitoepigenome as the nexus linking environmental exposures, mitochondrial (dys)function and altered transcriptional dynamics.

Given the importance of retrograde and anterograde communication between mitochondrial and nuclear genomes, mitoepigenomic modification remains an unexplored mechanism by which stress can be embedded at the transcriptional level via the mitochondrial response [4]. Certain environmental stressors and disease states are associated with alterations to the mitoepigenome [5,6,7] and have been used to investigate the extent and role of mitoepigenomic modifications, primarily in an in vitro context. Targeted bisulphite sequencing indicates that the mitochondrial D-loop control region is often differentially methylated upon exposure to, for example, airborne particulate matter [8]. An in vitro study using MeDIP-Seq to interrogate the entire mitochondrial genome, revealed that 7 mitochondrial genes are hypomethylated in hepatocytes exposed to valproic acid, including MT-CO1 and MT-CO2 [9]. At a whole organism level, there are a limited number of in vivo studies investigating mitoepigenomic response to environmental exposures. These studies often investigate a subset of mitochondrial genes, precluding a complete quantification of the mitoepigenetic response [8, 10, 11]. In addition to the paucity of data regarding chemical exposures and disease, it remains to be known whether non-chemical environmental stressors are able to affect the mitoepigenome. For example, it is currently unclear whether circadian disruption, a known affector of metabolism, influences the mitochondrial epigenome.

Circadian stress is an important model in which to investigate the mitoepigenomic response to stress. Circadian rhythmicity is essential in driving the appropriate temporal patterns of metabolism, meaning disruption of circadian rhythmicity results in aberrant metabolism [12]. Successful circadian rhythmicity requires a ‘zeitgeber’ to synchronise endogenous clocks, primarily in the form of photic input. In mammals the pineal gland is under the control of the suprachiasmatic nucleus (SCN) in the hypothalamus, which is affected by light exposures [13]. In birds, however, photic input occurs both through the retina and the pineal gland which functions as an independent oscillator and, together with the SCN, directs rhythmicity in behaviour, physiology and immunity in response to ambient light conditions [14, 15]. The pineal gland’s function lies primarily in the circadian secretion of melatonin, whose circulating levels vary in response to photoperiod. Deviation of light conditions from those required to generate a normal daily rhythm may disrupt the core circadian machinery, including the melatonin rhythmicity, leading to suboptimal physiological states [16,17,18]. Interestingly, melatonin is intimately linked to mitochondrial function, whose synthesis within the mitochondria may have arisen early in the endosymbiotic origins of eukaryotes [19]. As such, melatonin is thought to play a neuroprotective role, with circulating levels representative of age and susceptibility to psychiatric disorders. Indeed, melatonin’s protective action acting through mitochondria may provide an ameliorating role in the development of cancer cells [20]. Importantly, circadian stress is associated with alterations to methylation of the genomic DNA [21]. Whether such methylation changes are mirrored in the cell’s mitochondrial genome is yet to be determined.

Here we show that random illumination patterns elicit mito-epigenomic effects within the pineal gland and identify epigenetically plastic regions (mtDMRs) which respond to stress in a sex-specific manner. We make use of a MeDIP-Seq approach and a subsequent robust bioinformatic analysis, previously published by an independent group [22], that controls for the false detection of mitochondrial methylation resulting from nuclear mitochondrial pseudogenes. Using this method, we reveal the extent of mitochondrial cytosine methylation throughout the mitochondrial genome in the pineal glands of both male and female chickens (Gallus gallus domesticus) in response to externally induced circadian stress, taking advantage of the fact that avian pineal glands respond directly to photic stimuli.

Results

Sequencing and alignment

Sequencing of the 12 MeDIP-seq libraries relating to control and circadian stressed male and females gave an average of approximately 2.9 × 105 of reads assigned to the reference chromosome MT (Gallus gallus 6.0, NCBI). Upon controlling for nuclear mitochondrial pseudogenes, we mapped 2.7 × 105 reads, which corresponds to 2.2 ± 0.5 × 104 reads per individual, this corresponds to 93.6% of the reads mapped to the chrMT prior to correction. The average depth of these aligned reads was 172.6 X ± 39.2 X and covering 99.9% (± 0.1) of the chicken mtDNA genome. On average, the alignment against the mtDNA presented 32.6 times more depth of coverage than the whole chicken reference genome (Table S1).

Identification of differentially methylated CpNs and clustering within mtDMRs

Four pairwise comparisons were performed on the mtDNA sequencing data to identify differentially methylated CpN dinucleotides (DMCpNs) resulting from circadian stress and as a result of sex: Male control vs. female control (MC vs. FC), male stress vs. female stress (MS vs. FS), male control vs. male stress (MC vs. MS) and female control vs. female stress (FC vs. FS) (Fig. 1). We then clustered significant DMCpNs (P < 0.05) observed as peaks in Fig. 1 to define differentially methylated regions within the mitochondrial genome (mtDMRs) for downstream analysis.

Fig. 1
figure 1

Locations of DMCpNs within the mitochondrial genome. The solar plot displays the results of our differential methylation analysis, showing the significance of differential methylation for CpN dinucleotides across the mitochondrial genome for each of our pairwise statistical comparisons; Male control vs. female control (MC vs. FC), male stress vs. female stress (MS vs. FS), male control vs. male stress (MC vs. MS) and female control vs. female stress (FC vs. FS). Each data point represents a possible differentially methylated CpN dinucleotide and is positioned based on its base pair (BP) location on the mitochondrial genome and its significance (-log(p-value)). Dashed lines indicate p-values of 0.05 and 0.005. The shaded blue peak area represents an example of an mtDMR within the COX1 gene, in which CpN methylation is clustered

Full size image

We identified 62 mtDMRs across all statistical comparisons (Table S2), 46 of which had a length greater than 5 bp (range = 5-387 bp). 15 of these mtDMRs were present in at least two of the pairwise statistical comparisons, while 47 were unique regions within each contrast (Fig. 2A and C). Base composition within the mtDMRs did not differ significantly from the expected proportions within the mitochondrial genome, with the exception of G, which differed from the expected frequency in our MS vs. FS comparison (p = 0.01, Table S3). Methyl cytosines within the mtDMRs showed differential coverage in all four possible CpN dinucleotide combinations (table S4). The number of CpNs exhibiting differential coverage was not significantly different from the total number of each CpNs within a particular mtDMR region (Table S4). Furthermore, the proportions of CpNs within mtDMRs were not different from the expected proportions within the whole mitochondrial genome (Table S4). Although suggestive of differential cytosine methylation, differential coverage in MeDIP-Seq does not allow for the definitive identification of cytosine methylation with single base accuracy. It is therefore not possible to identify whether certain CpNs are more commonly differentially methylated than others based on differential coverage alone. With this in mind, an alternative way to determine the importance of different CpNs within the mtDNA is to look at the frequency of mtDMRs containing no instances of a particular CpN. In this case the absent CpN cannot be driving the significant differences in coverage, meaning that the differential coverage in these mtDMRs is a result of differential methylation of the other three CpNs. In this context, we found that CpG motifs were absent in 43.5% of the mtDMRs, whereas the other CpNs showed lower absences (CpA = 22.6, CpT = 17.7, CpC = 21).

Fig. 2
figure 2

The numbers and locations of regions of differential methylation (mtDMRs) in response to sex and stress. A The location of regions of differential cytosine methylation (mtDMRs) across the mitochondrial genome of pineal gland mitochondria resulting from circadian disruption (Stress). Within the 4 pairwise statistical comparisons, regions of hypo (blue) and hyper (red) methylation relate to our pairwise statistical comparisons; Male control vs. female control (MC vs. FC), male stress vs. female stress (MS vs. FS), male control vs. male stress (MC vs. MS) and female control vs. female stress (FC vs. FS). Coding regions within the heavy strand (mtDNA genes, above line) and light strand (mtDNA genes, below line) of the mitochondrial DNA are displayed for reference, with subunits of the same electron transport complexes shown in the same colour. tRNAs are shown in purple. Boxes around regions of differential methylation indicate the presence of sequence motifs associated with HNF4 (solid box), ATF4/ZNF324 (dashed box) or both (dot and dashed box). Letters within each window represent the nature of the mtDMR. a = sex-specific mtDMRs unrelated to stress, b = mtDMRs with identical responses in males and females, c = mtDMRs differing in the control and stress of one sex only, d = mtDMRs different across all statistical comparisons. B Venn diagram showing the mtDMRs which are lost, gained or remain between MC vs. FC and MS vs. FS, representing sex-specific stress-related regions (gained/lost) and sex-specific, stress-independent regions (remained). C A Venn diagram showing the numbers of mtDMRs which are either unique or present across our pairwise statistical comparisons

Full size image

Sex differences in mtDMRs

A comparison of mtDMRs present in the MC vs. FC with those present in MS vs. FS revealed a predominant increase in the number of mtDMRs upon exposure to circadian stress. A total of 9 mtDMRs (p-value ≤ 0.05), which were present in MC vs. FC comparison, were lost upon exposure to circadian stress (i.e. no longer present in the MS vs. FS) whilst 13 were gained, representing sex-specific stress-related regions (gained/lost). In comparison, five mtDMR sites remained, representing sex-specific, stress-independent regions (Fig. 2A and B). Of the total mtDMRs discovered, 25% were represented in at least two of the four contrasts. The extent of reoccurrence of mtDMRs in different contrasts allowed for the identification of the nature of their response to both stress and sex (Fig. 2A and C). We found that the most commonly recurring mtDMRs were those identified in both MS vs. FS and MS vs. MC comparisons (i.e. ‘Male stress-specific regions’, of which there were five mtDMRs), those differentially methylated only as a result of sex (i.e. ‘sex-specific regions’, which we found in both MC vs. FC and MS vs. FS comparisons, of which there were three mtDMRs) and regions which were differentially methylated as the result of sex and stressor (‘broadly-affected regions’, which were present in all statistical comparisons, of which there were two mtDMRs).

Hypo vs. hypermethylation

Within the mtDMRs, the direction of methylation was strongly influenced by stress and sex. Overall, mtDMRs were hypomethylated in FC in comparison to FS birds (10:4 hypo-:hyper- methylated mtDMRs in FC compared to FS) whereas MC vs. MS show equal levels of hyper and hypomethylation (Fig. 2A). The mtDMRs showed similar levels of hypo- and hyper-methylation when comparing MC vs. FC birds whereas MS birds showed a larger proportion of hypo- methylated regions compared to FS individuals (11:7 hypo-:hyper- methylated mtDMRs in MS compared to FS birds).

Motif discovery and enrichment analysis

Within mtDMRs, five significant motifs were found within the mtDMRs (XSTREME, E < 0.05, minimum width 6 bp, Table 1, Fig. 2A). Through a comparison of these motifs with three vertebrate motif databases, two of the motifs were found to share composition with known vertebrate motifs: The SGACTSTG motif was associated with the binding of hepatocyte nuclear factor alpha and gamma (HNF4A/HNF4G) and the CACCATCCT motif was associated with the binding of both zinc finger protein 324 (ZNF324) and the activating transcription factor 4 (ATF4). HNF4A/HNF4G binding-domains were identified in four mtDMRs located in all but the MC vs. FC comparison. Three of these were present in a single region from bases 6531–6625 (coding for tRNA-Cys and tRNA-Tyr) and were differentially methylated in MC vs. MS, FC vs. FS and MS vs. FS comparisons. Within this region, the pattern of methylation after stress differed in direction between males and females, being hypomethylated in stressed males (MS) in comparison to control males (MC) and hypermethylated in stressed females (FS) in comparison to their control counterparts (FC). Additionally, twelve ZNF324/ATF4 motifs were identified within mtDMRs and were found across all comparisons. Of these twelve sites, nine occurred within just 4 mtDMRs located within the ND2, COX1, ATP6 and CYTB coding regions. The broadly affected COX1 region (bp 7499–7717) showed the same pattern of hypomethylation in stressed males and females compared to their relative controls, although basal and stressed levels of methylation in females were higher than those of equivalent males. Other ATF4 sites commonly showed a pattern in which basal levels of methylation were unaffected by stress but were hypomethylated in females in comparison to males.

Table 1 Consensus motifs identified within mtDMRs of the mtDNA of chickens exposed to circadian stress. All significance values for the Homologous known vertebrate motifs identified through a comparison with vertebrate motif databases are shown where applicable. Sea P values were generated using XSTREME (Motif Discovery and Enrichment Analysis: Version 5.4.1 released on Sat Aug 21 19:23:23 2021 -0700) and represent the significance in the similarity between the sequence extracts and the motif sequence. The E value describes the same statistic accounting for multiple testing based on the number of identified motifs. It represents an estimate of the number of motifs that would be found with equal enrichment in shuffled versions of the positive sequences
Full size table

Motif presence across diverse taxa

Given the highly conserved nature of the mitochondrial genome across eukaryotes, we might expect that functionally important motifs within the mtDNA are conserved across diverse taxa. Having identified DNA-binding motifs (ATF4, ZNF324 and HNF4A) enriched in mtDMRs, we therefore sought to establish the prominence of these motifs within mitochondrial genomes of chickens (Gallus gallus), humans (Homo sapiens), alligators (Alligator mississippiensis), zebrafish (Danio rerio) and fruit flies (Drosophila melanogaster). Binding motifs were taken from human databases (tables S5, S6 and S7) and compared against the mitochondrial genomes. Comparisons were performed using FIMO, a web-based tool designed for this purpose by scanning sequences and determining sites which align with pre-defined motifs [23].

Identification of ATF4 recognition sites

Five instances of the ATF4 binding motif (Fig. 3A) were identified within the G. gallus genome, which were located in three protein-coding regions (16S rRNA, COX3 and ND4), while there were between two and four ATF4 sites present in the other species tested (H. sapiens (3), A. mississippiensis (2), D. rerio (4), and D. melanogaster (3)), with 76% of sites occurring on the mitochondrial light (negative) strand (Fig. 3B, Table S8). Two protein-coding regions, 16S ribosomal rRNA and ND4 contained 53% of these sites and one region was found within the 16S ribosomal RNA gene in all species except for D. melanogaster. Alignment of the common ATF4 binding domains within the 16 s RNA sequences (± 20 bp upstream and downstream) revealed 100% conservation between species in the core binding site sequence (TGTTGGATCAGGAC) and a 70% consensus sequence differing in only 2/54 bases (Fig. 3C).

Fig. 3
figure 3

ATF4 presence in mtDNA across diverse taxa. A The DNA recognition motif for human ATF4 transcription factor used for a comparison across diverse taxa. B Locations of human ATF4 transcription factor recognition sites identified using FIMO (black bars above mtDNA coding regions) within the mitochondrial genomes of diverse taxa: Gallus gallus, Homo sapiens, Alligator mississippiensis, Danio rerio and Drosophila melanogaster. Coding regions within the heavy strand (mtDNA genes, above line) and light strand (mtDNA genes, below line) of the mitochondrial DNA are displayed for reference, with subunits of the same electron transport complexes shown in the same colour. tRNAs are shown in purple. C Sequence allignment of the 16S rRNA region associated with ATF4 binding (shown with star in B, outlined in dashes in C) as well as 20 bp upstream and downstream of the site

Full size image

Identification of ZNF324 recognition sites

The chicken mitochondrial genome contained 10 sites sharing sequence with that of the ZNF324 recognition site (Figure S1, Table S9). Despite being found in all species tested, fewer sites were located in the mitochondrial genomes of H. sapiens (6), A. mississippiensis (8), D. rerio (5), and D. melanogaster (2). Recognition sites were found in the ND2, ND4 and ND5 genes in three of the five species, each coding for components of electron transport chain complex I. 94% of the recognition sites were located within the mitochondrial light (negative) strand.

Identification of HNF4A recognition sites

Our analysis revealed five regions of the G. gallus mtDNA containing sequences matching the HNF4A recognition domain (Figure S1, Table S10), with fewer sites found in H. sapiens (4), A. mississippiensis (2), D. rerio (1), and D. melanogaster (0). Unlike ZNF324 and ATF4 sites, location of the HNF4 recognition sites was variable across species with no apparent commonalities in location. 92% of sites were located on the mitochondrial heavy (positive) strand.

Discussion

With few exceptions, both wild and domesticated species have evolved under the influence of circadian fluctuations in light intensity. Disruption of the resultant evolved circadian rhythmicity represents a form of environmental stress that is relevant in multiple contexts from animal welfare to human health and ecology [24, 25]. Our results demonstrate that the mitoepigenome of the pineal gland, a core regulator of circadian rhythmicity in vertebrates, is impacted by circadian stress during postnatal development when deviations from a 12:12 light dark cycle occur early in life. Importantly, the nature of this epigenetic response differs between the sexes.

The effects of circadian stress upon mtDNA methylation

The notion of mtDNA methylation remains controversial, with numerous studies reporting an absence or minimal level of methylation [26,27,28,29]. Our data clearly demonstrate that in the pineal gland of chickens not only is mtDNA methylated, but cytosinic methylation of mtDNA is variable upon developmental exposure of animals to circadian stress. Moreover, the effects are different between the sexes. The novel finding that mtDNA within the pineal gland is extensively methylated, often in non-CpG contexts within distinct epigenetically plastic regions (mtDMRs) of the mitochondrial genome has stark implications for our understanding of mitochondrial function and malfunction. Although circadian disruption is known to have metabolic effects, this study provides the first evidence that an environmental stressor, as opposed to a disease state [2, 6] or chemical exposure [8,9,10,11], is capable of modifying the epigenetic landscape of the mitochondria. Additionally, the identification of transcription factor binding motifs within the mtDMRs identified adds evidence to the notion [30] that nuclear-derived transcription factors play a role in regulating mtDNA transcription across species.

Locations of mtDMRs

Although mtDMRs were found throughout the mitochondrial genome, across all statistical comparison, mtDMRs were biased towards specific coding regions of the mtDNA. For example, CYTB, ND1, COX1 and ND5 contained over 50% of all the detected mtDMRs (Fig. 2A). All four regions are protein coding genes for key components of electron transport chain respiratory complexes I (ND2, ND4), complex III (CYTB) and IV (COX1). Altered transcriptional activity of these genes would have clear functional consequences for electron transport, proton gradient formation and thus the ATP production efficiency of the affected mitochondria.

Overall, the finding that specific regions of the mtDNA are more responsive than others lends strong evidence to the notion that the mitoepigenome is regulated and responsive to the environment in a targeted manner. The observation that the CpN composition within mtDMRs is not different from that of the mitochondrial genome as a whole adds weight to this notion, suggesting specific motifs may be important in driving mtDMR methylation as opposed to cytosine content, as is the case in CpG islands within the nucleus [31].

Certain genes within the mtDNA contained very few to no mtDMRs, indicating a lack of cytosinic methylation. Of particular interest is the finding that cytosine methylation of the D-loop control region is not responsive to stress. The D-loop is central to mitochondrial transcription and mtDNA replication, representing the origin of transcription around the DNA plasmid. Because of its importance in transcriptional dynamics, previous studies have often focused on the D-loop in order to detect methylation. Alterations to D-loop methylation have been previously observed in response to environmental stress [8, 10, 11, 32] and disease [33,34,35,36]. The absence of D-loop mtDMRs here may indicate that this region is less responsive to circadian stress. Alternatively, it is possible that individual cytosines at specific locations within the D-loop are important but were missed in our analysis, which was biased towards regions of differential methylation spanning multiple cytosine bases.

Sex differences in mtDMRs

Despite the potential for the mitochondrial epigenome to influence mitochondrial function, mtDNA has been seldom investigated in relation to sex differences in metabolism. Our data reveal a clear sex difference in the mitoepigenome, both in terms of baseline methylation levels between control male and female birds and in the response of the pineal gland to circadian stress. By comparing the differential methylation of mtDMRs between control and stressed male and females, we were able to identify regions of the mtDNA associated with sex, stress and the interaction of the two. Control males and females differed in mtDNA methylation within one tRNA, 2 rRNA and 7 protein-coding genes, 3 of which coded for sub-units of electron transport chain complex 1 (ND1, ND2 and ND5). Interestingly mtDMRs within atp6, trna1 and nd2 remained differentially methylated between both stressed and unstressed males and females, representing sex-related epigenetic elements. Furthermore, two of these ‘sex-specific regions’ contained binding sites for the Atf4 transcription factor, hinting at a basal sex difference in mito-nuclear interactions. These findings are concomitant with findings showing sex differences in metabolism and mitochondrial function throughout the animal kingdom, with both animal and human studies demonstrating divergence in the stress responses of males and females [37,38,39]. Indeed, white leghorn chickens display clear sexual dimorphism and differ in their basal rates of metabolism as well as in their physiological stress responses [40].

Previous studies in pigs indicate opposite directionality of mtDNA promoter methylation and transcription between males and females in response to stress [32]. Looking at the entire mitochondrial genome, we find elements in which mtDNA which respond both oppositely (e.g. trna10, which is hypomethylated in stressed males and hypermethylated in stressed females) and similarly (cox1, which is hypomethylated in stressed males and females), yielding a more complex picture of the mitochondrial methylome. Interestingly, the trna10 site that responds oppositely contains a motif recognition site for Hnf4 stress-related transcription factor whereas the cox1 mtDMR contains an Atf4 stress-related transcription factor binding site.

Should these transcription factors bind differentially based on methylation state, these results hint at a complex picture in which different nuclear-derived transcription factors exert differential effects based both upon sex and the location within the mtDNA.

Cytosine context within mtDMRs

In relation to the nucleotides adjacent to methyl-cytosines within the mtDNA, our results agree with previous findings from diverse taxa such as mice [41] and fish [42], suggesting that non-CpG methylation is relevant and widespread in vertebrates. Investigations utilizing NGS with primarily mammalian models have detected higher non-CpG vs CpG methylation within the mtDNA (e.g. [41, 43]). Non-CpG methylation is more commonly observed in prokaryotes, however is also relevant to eukaryotic genomic DNA (despite the predominance of CpG methylation) and is increasingly being considered as relevant to transcriptional processes within the nucleus [43, 44]. As opposed to CpG islands found within the nuclear genome, significant regions of differential methylation within the mtDNA contained methylated cytosines within CpN motifs in similar proportions to those within the mtDNA as a whole. Hence, the location of the methylated region within specific mtDNA motifs may be more critical to any functionality or DNA response than simply CpN bases being methylated per se. Methylated cytosines are more likely to transition to a T during DNA replication of the DNA. As a result of the predominance of methyl cytosines in a CpG context within the nucleus, it is perhaps not surprising that across the animal kingdom, age-related C > T transitions within the genomic DNA occur primarily in CG dinucleotides [45]. Interestingly, the data of Cagan et al., 2022 reveal that C > T transitions in mtDNA with age occur commonly within a CpT context (particularly at GCT motifs), in keeping with our data showing extensive CpT representation within regions of differential methylation. mtDNA methylation may therefore not be only of interest in a mechanistic context but also in an evolutionary context, in which C > T transitions within certain regions may have shaped the mitochondrial genome through time [46]. It would therefore be of interest to determine whether C > T transitions within the mtDNA are located within distinct regions across genera.

Mechanistic consequences of differential methylation

Having discovered extensive cytosine methylation within the mtDNA, it is important to establish the potential mechanistic consequences of differential methylation within the mtDNA. By identifying functional regions which are plastic in their epigenetic response, we suggest that differential cytosine methylation within the mtDNA is indicative of functionality. MtDNA transcription is understood to start from the D-loop control region and the resultant polycistronic transcript is separated via excision of the tRNA prior to post-transcriptional processing [47]. Although mitochondrial gene promoters are thought to lie solely within the control region of the mtDNA, it has been suggested that a plethora of nuclear-derived DNA-binding factors are responsible for the regulation of mtDNA transcription elsewhere within the genome [30, 48, 49].

Within the nucleus, DNA methylation has been associated with the repression of transcription factor binding [50]. However, the regulation of methylation may be more complex, with suggestions that transcription factor binding may, itself, modulate DNA methylation state [51]. Regardless, the finding here that methylation varies across treatments and sexes in specific parts of the mitochondrial genome might indicate an important role of cytosine modification in the regulation of specific mitochondrial genes, a hypothesis that requires further investigation.

Here, mtDMRs responding to circadian stress contained binding motifs associated with nuclear-derived transcription factors known to be central to metabolic homeostasis and the cellular stress responses (HNF4, ATF4 and ZNF324). For example, HNF4A is important in regulating metabolism, exerting an influence over glycolysis and mitochondrial metabolism and has been demonstrated to bind directly to mtDNA [52]. The importance of such transcription factors in modulating stress-responses coupled with the increasing evidence that nuclear transcription factors locate to the mitochondria [30] add weight to the novel idea that nuclear stress signals directly alter mitochondrial gene transcription through binding of DNA elements throughout the mitochondrial genome. It appears that mtDNA transcription is therefore driven by more than simply the mitochondria-specific POLRMT polymerase and the TFAM and TFB2M transcription factors thought to dominate mtDNA transcription.

Among the transcription factors found to have binding motifs within mtDMRs, ATF4 represents a key component of the mitochondrial unfolded protein stress response and integrated stress response [53]. Importantly, although there have been relatively few studies investigating nuclear transcription factor binding to mtDNA, one of ATF4’s orthologues in C. elegans, ATFS-1 has been shown to bind mtDNA, modifying transcription of oxidative phoshporilation genes [54]. ATF4 is correlated with mitochondrial function in non-stressed cells, suggesting a role in not only the mitochondrial stress response, but in normal mitochondrial homeostasis [55]. For example, exposure of HeLa cells to 4 differing mitochondrial stressors revealed that ATF4 binding sites were present in many of the genes involved in the varying stress responses [55]. Crucially, mitochondrial stress was found to inhibit mitochondrial translation and decrease mitochondrial ribosomal proteins in what the authors state as a ‘localized response’. Our findings of a highly conserved region across species within the 16 s ribosomal DNA of the mitochondria lends evidence to the idea that ATF4 is, in fact, be mechanistically linked to the alterations in translation observed in stressed mitochondria as opposed to part of a separate response.

As well as identifying motifs consistent with DNA-binding factors, our mtDMRs reveal additional novel motifs with no clear eukaryotic homologues. Given that, unlike the nuclear DNA, mtDNA is bacterial in origin, it should also be considered that mtDNA transcription could additionally be regulated in a ‘non-nuclear’ way. Bacterial plasmids are actively methylated and regulatory elements such as enhancers are found in bacteria (Xu and hoover, 2001). The mtDMRs detected without eukaryotic equivalents may therefore represent important motifs involved in an, as yet, undetermined process of mtDNA regulation.

Understandably, the close link between transcription and methylation within the nucleus bias functional explanations of mtDNA methylation towards transcription factor binding. Despite this, there are alternative functional explanations. DNA methylation, through the association with methylated DNA-binding proteins, can influence the rate of polymerase translocation across the DNA. mtDMRs could thus represent pause or slow sites across the mtDNA akin to those important for gene regulation in bacteria [56]. Methylation of the mtDNA could also be involved in altering DNA stability and 3-dimensional structure. In the absence of histone proteins, mtDNA accessibility can be altered by the 3-dimensional structure of the nucleoid, altering transcriptional dynamics. It is realistic to hypothesize that mtDNA methylation could influence the conformation of the mtDNA or, by altering the affinity of DNA binding proteins such as TFAM, could influence DNA packaging and flexibility [57, 58].

Taking the present data together, it is possible to make predictions regarding the transcriptional responses of pineal cell mitochondria to circadian stress under a scenario in which ATF4, ZNF324 and HNF4 could bind specific regions of the mtDNA depending on their cytosine methylation status. The methylation-dependent binding of these transcription factors could modify transcription. For example the ND2, ATP6 and CYTB genes contain sites which are hypermethylated in male birds in comparison to female birds, regardless of whether they were subjected to the experimental stress or not. We might therefore expect parallel differences in transcriptional patterns when comparing these genes across the sexes, for example with relative repression of transcription in males should methylation repress transcription as in the nuclear genome. Conversely, COX1 may show opposite transcriptional responses across the sexes as it contains regions which are hypomethylated in male birds in comparison to female birds. This same region is hypermethylated in control birds in comparison to the sex-matched stressed birds, suggesting a similar transcriptional response in both stressed males and stressed females. Given the centrality of all of the protein coding genes to the dynamics of mitochondrial electron transport, methylation fingerprints within the mtDNA might be associated with alteration of transcription factor binding. Methylation-dependent altered binding of these transcription factors could result in changes in mitochondrial function via the transcriptional regulation of specific mitochondrial genes potentially regulated by these transcription factors.

Evolutionary perspective

Transcription factor binding specificities are largely conserved across the bilateria [44, 52, 59] and the roles of ATF4, HNF4 and zinc finger proteins as regulators of metabolism, as well as their expression patterns, are highly conserved from flies to mammals [60]. Therefore, given our observation of recognition sites for these DNA-binding proteins across a range of metazoan species, it is reasonable to infer that mitochondrial roles of these nuclear transcription factors may also have evolved early in metazoans. In addition, the identification of binding sites through mtDMRs suggest that methylation may be an important regulator of this transcription factor binding, a known phenomenon regarding protein-DNA interactions [61, 62] for example involving zinc-finger proteins [63]. Importantly, the presence of extensive non-CpG methylation within genomic DNA [64, 65] together with the observation that CpA methylation is important in altering DNA binding affinity [62, 63] suggest an important role of non-CpG methylation within the nucleus as well as within the mitochondria.

Conclusions

In conclusion, our results demonstrate that the mitoepigenome within pineal cells of chickens is responsive to circadian stress and varies with sex, adding weight to evidence showing that mtDNA methylation could be one the one hand responsive to environmental conditions and on the other hand important for the transcriptional regulation of specific mitochondrial genes independent of the D-loop. We have identified regions of differential methylation within the mtDNA that appear to have a role in the nuclear-driven regulation of mitochondrial transcription in response to environmental stress. Although it is not possible to determine to what extent mtDNA methylation represents a causative or correlative element, these data suggests that the binding of the nuclear factors ATF4 and HNF4 to the mtDNA could be relevant for mitochondrial function and dependent on mt DNA methylation levels. In order to unpick this response, it is essential that future experiments aim to connect mtDNA methylation, transcription factor binding and differential transcription of the mitochondrial and nuclear genomes.

Methods

Animals and husbandry

Fertilized Hy-Line (Swedfarm, Linghem, Sweden) white leghorn chicken (Gallus gallus domesticus) eggs were incubated in a 60% humidified atmosphere at 37.8 °C with rotation once every hour, until hatching. After hatching, all individuals were wing tagged, placed in boxes, and directly transferred to cages. They lived in mixed groups containing both sexes (females n = 18, males n = 16), and between 5 and 6 chickens per cage. They were not confidently sexed until culling and were therefore not balanced between groups. Throughout the experiment all birds had ad libitum access to water and food as well as wood chips on the floor of every cage. Feed trays including starter food provided access to food within the first days, which was later on replaced by hanging food (Pennafood). During the first period of the experiment, chickens had access to heated roofs. The animals were raised in temperature and light-emitting diode (LED) light-controlled cages.

Treatment

The animals were randomly divided into two groups after one day of incubation. The control group (n = 6, 3 males and 3 females) was kept at standard 12:12 light–dark cycle for their entire life. The circadian stress group (n = 6, 3 males and 3 females) was exposed to continuous light from days 3–6 of age. From days 7–24, birds experienced a randomized light schedule (36 h of light every 3 days delivered at random intervals between 3 to 21 h). Over a three-day period, the total amount of light and dark was equal to that of control animals. From day 25 and for the rest of the experiment, the animals were kept at standard 12:12 light–dark cycle.

Each cage received light from two LED lamps, controlled by a timer. Cages were isolated to stop light from reaching neighbouring cages, confirmed by light measurements before the experiments. Culling took place over a three-day period when the birds were 34–36 days old. Control and circadian-stressed chicks were sacrificed by decapitation and pineal glands were removed from the brain and snap frozen in liquid nitrogen. Data from all individuals were included in the final analyses without any exclusions.

DNA isolation

In order to isolate total DNA from pineal gland tissue, Allprep® DNA/RNA/miRNA Universal kit (Qiagen GmbH, Hilden, Germany, Catalog # 80,224) was used. First, the complete pineal gland was diluted in 350 μl RLT Plus buffer (Qiagen GmbH, Hilden, Germany, Catalog # 1,053,393), which contained guanidine thiocyanate for tissue lysis and then homogenized for 40 s in a rotor–stator homogenizer (FastPrep−24 Homogenizer, MP™ Biomedicals, California, USA). Further steps were performed according to the manufactures protocol. Samples was stored at -80 °C. The DNA concentration was determined by nanodrop spectrophotometer and calculated by Gene Quant pro RNA/DNA calculator Software (NanoDrop ND−1000 Spectrophotometer, Thermo Scientific, Wilmington, USA).

Ion torrent MeDIP-Seq library preparation

MeDIP-Seq libraries for Ion Torrent semiconductor sequencing were performed on the Ion Torrent PGM using a modified protocol from the Ion Plus Fragment Library kit (ThermoFisher, Massachusetts, USA, Catalog # 4,471,252) according to a previously optimized protocol (Guerrero-Bosagna & Jensen, 2015). For 5mC MeDIP-Seq, 8 μg of DNA was sheared to a mean fragment size of 350 bp by sonication in a Fisher ultra-sonicator, with probe attached to a cooling cup horn. The final preparation of the libraries was performed at the SNP&SEQ facilities of the SciLifeLab (Sweden) and consisted of the end-repair, purification with Agencourt AMPure XP reagent beads (Beckman Coulter, California, USA, Product # A63880), and ligation with Ion Torrent compatible sequencing adapters following the standard Ion Plus Fragment Library Kit protocol. Adapter-ligated DNA was purified using Agencourt AMPure XP reagent beads and immunoprecipitated using a 5-methylcytosine kit (MagMeDIP, Diagenode, Liege, Belgium, Catalog # C02010021) according to the manufacturer's protocol. Following 5mC MeDIP, libraries were size-selected from 200–400 bp on a 2% agarose E-gel (Invitrogen, Massachusetts, USA, Catalog # A45205) and extracted using a MinElute Gel DNA extraction kit (Qiagen GmbH, Hilden, Germany, Catalog # 28,604) according to the manufacturer's instructions. Libraries were then amplified at 18 cycles of PCR (95 °C, 15 s; 58 °C, 15 s; 70 °C, 1 min), purified using MinElute columns (Qiagen, Hilden, Germany, Catalog # 28,004), and eluted in 16 μl elution buffer (Qiagen, Hilden, Germany, Catalog # 19,086). Libraries were assessed for quality and concentration on an Agilent Bioanalyzer instrument.

Ion torrent sequencing

MeDIP libraries were templated using the Ion OneTouch 2 System instrument with the Ion PI Template OT2 Kit v2 for the Ion Proton instrument according to manufacturer's protocols (Life Technologies). Following template reactions, templated libraries were assessed for quality using the Ion Sphere Quality Control Kit (Catalog # 4,468,656) and Qubit fluorometer (Life Technologies). Templated libraries that passed manufacturer recommended criteria were enriched using the Ion OneTouch ES instrument. Immediately after enrichment, MeDIP libraries were sequenced on the Ion Proton with PI chips using the Ion PGM 200 Sequencing kit or Ion PtdIns Sequencing 200 Kit v2, according to instructions provided by the manufacturer (Life Technologies).

MeDIP-Seq analysis

Ion Torrent sequencing reads were quality-filtered and aligned to the chicken reference sequence Gallus_gallus 6.0 (NCBI) using Torrent Suite Software's (v4.0) TMAP alignment software. Alignment files were processed into the BAM format and used for further analyses (Figure S2). The index and coverage depth of the “.bam” files were performed using Samtools v.1.11 [66] with the “index and “depth” options.

To control for nuclear-mitochondrial pseudogenes within Methylated DNA Immunoprecipitation Sequencing data, we used a recently published pipeline [22], which consists in four steps using the previously mentioned BAM files as input. First, we extracted the new BAM file of “unique reads”, which were aligned to no more than one specific region of the reference genome (RG). For that, we used Samtools v.1.11. Second, we extracted a FASTQ file with those “unique reads” from the previous generated BAM files using BEDTools v.2.29.2 to be aligned against an RG index without the chrMT. This reference genome was build using “bowtie2-build” function from Bowtie v.2–2.4.2. Third, we extracted a new BAM file containing the non-aligned reads from the previous BAM file generated by the last step. Forth, we extract the.fasta file from the bam file of the earlier step to be aligned to a new genome index, but now build using only the mtDNA sequence from the RG.

Through these last BAM files generated from the previous pipeline, the MEDIPS R (Lienhard et al., 2014) package was used for basic processing, quality controls, normalization, to assess the proportion of genome-wide CpGs covered in each sample and identification of differential coverage between different contrasts. We used the package's default parameters, as we have previously described (Pértille et al., 2017), except for the windows size parameter in which we mapped the CpN positions of the chicken reference genome to be used as region of interest (ROI), and also the chromosome selection parameter, of which we used only the chicken mtDNA. To map the CpNs we used the matchPattern function from “stringr” package from R. BSgenome.Ggallus.UCSC.galGal6 package from Bioconductor repository was uploaded as the reference of the chicken mitochondrial genome. After that, we create a function in R to merge, every single CpN tested for differential methylation with its respective mt annotation coordinates obtained from the Variant Effect Predictor (VEP) tool v.71 (McLaren et al., 2010). This function used part of the script available at https://github.com/stephenturner/solarplot/blob/master/solarplot.R. Moreover we also used plot functions in the R environment to show the differential methylation statistics distributed across the contrasts we have tested.

Differential methylation analysis

For the 12 individuals, an average “enrichment score” of 2.2 (± 0.1) was obtained to test the CpGs enrichment within the genomic regions covered by the given set of short reads compared to the full reference genome for which the “enrichment score” was 1.1 (table S11). This score was obtained dividing the relative frequency of CpGs within the given regions by the reference genome sequences. Non-enriched DNA sequences should present values close to 1, in contrast, a MeDIP-Seq experiment should return CpG rich sequences that will be indicated by increased CpG enrichment scores (Lienhard et al., 2014). We then performed a series of differential methylation analyses for each CpN combination in order to identify differentially methylated CpNs (DMCpNs) between our four statistical contrasts. Our method interrogated differential coverage within each contrast, individually, for each one of the CpN sites present on the chicken mtDNA genome that is 16,775 bps wide. For each site, we tested four contrasts: two comparing control vs stressed groups for males and females separately (MC vs MS and FC vs FS, respectively) and two to test for sex differences between the control and stressed groups (MC vs FC and MS vs FS, Fig. 1). The analysis of differential methylation, using the MEDIPS package in R, utilizes a weighted trimmed mean of the logarithm of expression ratios (TMM) [67, 68]. By default, a minimum sum of counts across all samples per window of 10 is employed. MEDIPS uses the edgeR statistical method and we used a p-value threshold of ≤ 0.05 to define a significant DMR [67].

Identifying mtDMRs

Upon visualising the significant DMCpNs obtained through our contrasts, we found that DMCpNs clustered in distinct regions within the mtDNA, Fig. 1, supplementary table S2). These clusters of DMCpNs were delimited for downstream analysis, and defined as mtDMRs. For this clustering we constructed a model in which an mtDMR was defined by having a maximum gap of length x base pairs between adjacent DMCpNs. Through checking the number of merged windows generated with gap lengths, x ranging from 1–1962, we identified four important x-value ranges (127–179, 210–284,389–522 and 659–882). Within each range, increasing the gap length did not alter the number of merged regions. From 127 to 179, we identified the same 21 regions, giving a consecutive index (CI = 53, representing the range of x for which the number of significant regions remains constant), while from 210 to 284, we identified 16 merged windows (CI = 76), from 389 to 522 we identified 11 windows (CI = 124), and finally from 659 to 882, we identified 6 windows (CI = 223). High CI indicates low specificity because large changes in gap length do not change the number of dinucleotides which belong to that specific region. On the other hand, a low CI tends towards the border of significance among dinucleotides. We therefore used the smallest gap length (x = 127) that gave us the lowest CI, giving a conservative estimate of the number of mtDMRs.

Motif discovery and enrichment analysis

In order to identify motifs within mtDMRs as well as elucidate their functionality, 46 of the of the methylation-enriched regions (5–387 bp in length) were submitted to XSTREME motif discovery and enrichment analysis, which both detects enriched motifs and then searches for homologues within vertebrate databases (version 5.5.0 [69], E values ≤ 0.05). Having identified enriched motifs within our mtDMRs and any vertebrate homologues within them, the sequence of homologous motifs were then drawn from human databases (for motifs, see tables S6, S7 and S8) and the occurrences of these sequences were determined within the mitochondrial genomes of Gallus gallus, Homo sapiens, Alligator mississippiensis and Danio rerio using FIMO (version 5.5.0 [23], p < 0.0001).

Availability of data and materials

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The datasets generated for this study can be found in the ENA repository project PRJEB59721 (https://www.ebi.ac.uk/ena/browser/view/PRJEB35831).

References

  1. van der Wijst MGP, Rots MG. Mitochondrial epigenetics: an overlooked layer of regulation? Trends Genet. 2015;31(7):353–6.

    PubMed  Google Scholar 

  2. Devall M, Roubroeks J, Mill J, Weedon M, Lunnon K. Epigenetic regulation of mitochondrial function in neurodegenerative disease: New insights from advances in genomic technologies. Neurosci Lett. 2016;625:47–55.

    CAS  PubMed  PubMed Central  Google Scholar 

  3. Dimauro S, Davidzon G. Mitochondrial DNA and disease. Ann Med. 2005;37(3):222–32.

    CAS  PubMed  Google Scholar 

  4. Picard M, Turnbull DM. Linking the Metabolic State and Mitochondrial DNA in Chronic Disease, Health, and Aging. Diabetes. 2013;62(3):672–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  5. Lambertini L, Byun HM. Mitochondrial epigenetics and environmental exposure. Curr Environ Health Rep. 2016;3(3):214–24.

    PubMed  Google Scholar 

  6. Stoccoro A, Coppedè F. Mitochondrial dna methylation and human diseases. Int J Mol Sci. 2021;22(9):1–27.

    Google Scholar 

  7. Iacobazzi V, Castegna A, Infantino V, Andria G. Mitochondrial DNA methylation as a next-generation biomarker and diagnostic tool. Mol Genet Metab. 2013;110(1–2):25–34.

    CAS  PubMed  Google Scholar 

  8. Janssen BG, Byun HM, Gyselaers W, Lefebvre W, Baccarelli AA, Nawrot TS. Placental mitochondrial methylation and exposure to airborne particulate matter in the early life environment: An ENVIRONAGE birth cohort study. Epigenetics. 2015;10(6):536–44.

    PubMed  PubMed Central  Google Scholar 

  9. Wolters JEJ, Van Breda SGJ, Caiment F, Claessen SM, De Kok TMCM, Kleinjans JCS. Nuclear and Mitochondrial DNA Methylation Patterns Induced by Valproic Acid in Human Hepatocytes. Chem Res Toxicol. 2017;30(10):1847–54.

    CAS  PubMed  PubMed Central  Google Scholar 

  10. Armstrong DA, Green BB, Blair BA, Guerin DJ, Litzky JF, Chavan NR, et al. Maternal smoking during pregnancy is associated with mitochondrial DNA methylation. Environ Epigen. 2016;2(3):dvw020-dvw.

    Google Scholar 

  11. Liu Q, Li H, Guo L, Chen Q, Gao X, Li Ph, et al. Effects of short-term personal exposure to air pollution on platelet mitochondrial DNA methylation levels and the potential mitigation by L-arginine supplementation. J Hazard Mater. 2021;417(January):125963.

    CAS  Google Scholar 

  12. Fatima N, Rana S. Metabolic implications of circadian disruption. Pflueg Arch Eur J Physiol. 2020;472(5):513–26.

    CAS  Google Scholar 

  13. Ko CH, Takahashi JS. Molecular components of the mammalian circadian clock. Hum Mol Genet. 2006;15(2):R271–7.

  14. Cassone VM, Kumar V. Circadian Rhythms. In: Sturkie’s avian physiology (sixth edition). Scanes: Academic Press; 2015.

  15. Karaganis SP, Bartell PA, Shende VR, Moore AF, Cassone VM. Modulation of metabolic and clock gene mRNA rhythms by pineal and retinal circadian oscillators. Gen Comp Endocrinol. 2009;161(2):179–92.

    CAS  PubMed  Google Scholar 

  16. Gooley JJ, Chamberlain K, Smith KA, Khalsa SB, Rajaratnam SM, Van Reen E, et al. Exposure to room light before bedtime suppresses melatonin onset and shortens melatonin duration in humans. J Clin Endocrinol Metab. 2011;96(3):E463–72.

    CAS  PubMed  Google Scholar 

  17. Moaraf S, Vistoropsky Y, Pozner T, Heiblum R, Okuliarová M, Zeman M, et al. Artificial light at night affects brain plasticity and melatonin in birds. Neurosci Lett. 2020;716:134639.

    CAS  Google Scholar 

  18. Plano SA, Casiraghi LP, García Moro P, Paladino N, Golombek DA, Chiesa JJ. Circadian and metabolic effects of light: implications in weight homeostasis and health. Front Neurol. 2017;8:558.

    PubMed  PubMed Central  Google Scholar 

  19. Tan D-X, Manchester LC, Liu X, Rosales-Corral SA, Acuna-Castroviejo D, Reiter RJ. Mitochondria and chloroplasts as the original sites of melatonin synthesis: a hypothesis related to melatonin’s primary function and evolution in eukaryotes. J Pineal Res. 2013;54(2):127–38.

    CAS  PubMed  Google Scholar 

  20. de Almeida Chuffa LG, Seiva FRF, Cucielo MS, Silveira HS, Reiter RJ, Lupi LA. Mitochondrial functions and melatonin: a tour of the reproductive cancers. Cell Mol Life Sci. 2019;76(5):837–63.

    PubMed  Google Scholar 

  21. Cedernaes J, Schönke M, Westholm JO, Mi J, Chibalin A, Voisin S, et al. Acute sleep loss results in tissue-specific alterations in genome-wide DNA methylation state and metabolic fuel utilization in humans. Sci Adv. 2018;4(8):eaar8590.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. Devall M, Smith RG, Jeffries A, Hannon E, Davies MN, Schalkwyk L, et al. Regional differences in mitochondrial DNA methylation in human post-mortem brain tissue. Clin Epigenet. 2017;9(1):1–15.

    Google Scholar 

  23. Grant CE, Bailey TL, Noble WS. FIMO: scanning for occurrences of a given motif. Bioinformatics. 2011;27(7):1017–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. Dominoni DM, Borniger JC, Nelson RJ. Light at night, clocks and health: from humans to wild organisms. Biol Lett. 2016;12(2):20160015.

    PubMed  PubMed Central  Google Scholar 

  25. Bittman EL, Kilduff TS, Kriegsfeld LJ, Szymusiak R, Toth LA, Turek FW. Animal care practices in experiments on biological rhythms and sleep: report of the Joint Task Force of the Society for Research on Biological Rhythms and the Sleep Research Society. J Am Assoc Lab Anim Sci. 2013;52(4):437–43.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Mechta M, Ingerslev LR, Fabre O, Picard M, Barrès R. Evidence suggesting absence of mitochondrial DNA methylation. Front Gen. 2017;8(NOV):1–9.

    Google Scholar 

  27. Hong EE, Okitsu CY, Smith AD, Hsieh C-L. Regionally specific and genome-wide analyses conclusively demonstrate the absence of CpG methylation in human mitochondrial DNA. Mol Cell Biol. 2013;33(14):2683–90.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. Matsuda S, Yasukawa T, Sakaguchi Y, Ichiyanagi K, Unoki M, Gotoh K, et al. Accurate estimation of 5-methylcytosine in mammalian mitochondrial DNA. Sci Rep. 2018;8(1):5801.

    PubMed  PubMed Central  Google Scholar 

  29. Guitton R, Nido GS, Tzoulis C. No evidence of extensive non-CpG methylation in mtDNA. Nucleic Acids Res. 2022;50(16):9190–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. Marinov GK, Wang YE, Chan D, Wold BJ. Evidence for Site-Specific Occupancy of the Mitochondrial Genome by Nuclear Transcription Factors. PLoS ONE. 2014;9(1): e84713.

    PubMed  PubMed Central  Google Scholar 

  31. Gardiner-Garden M, Frommer M. CpG islands in vertebrate genomes. J Mol Biol. 1987;196(2):261–82.

    CAS  PubMed  Google Scholar 

  32. Jia Y, Li R, Cong R, Yang X, Sun Q, Parvizi N, et al. Maternal low-protein diet affects epigenetic regulation of hepatic mitochondrial DNA transcription in a sex-specific manner in newborn piglets associated with gr binding to its promoter. PLoS One. 2013;8(5):e63855.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. Stoccoro A, Siciliano G, Migliore L, Coppedè F. Decreased Methylation of the Mitochondrial D-Loop Region in Late-Onset Alzheimer’s Disease. J Alzheimer’s Dis. 2017;59(2):559–64.

    CAS  Google Scholar 

  34. Tong H, Zhang L, Gao J, Wen S, Zhou H, Feng S. Methylation of mitochondrial DNA displacement loop region regulates mitochondrial copy number in colorectal cancer. Mol Med Rep. 2017;16(4):5347–53.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. Zheng LD, Linarelli LE, Liu L, Wall SS, Greenawald MH, Seidel RW, et al. Insulin resistance is associated with epigenetic and genetic regulation of mitochondrial DNA in obese humans. Clin Epigenet. 2015;7(1):1–9.

    CAS  Google Scholar 

  36. Mishra M, Kowluru RA. Epigenetic modification of mitochondrial DNA in the development of diabetic retinopathy. Investig Ophthalmol Vis Sci. 2015;56(9):5133–42.

    CAS  Google Scholar 

  37. Demarest TG, McCarthy MM. Sex differences in mitochondrial (dys)function: Implications for neuroprotection. J Bioenerg Biomembr. 2015;47(1–2):173–88.

    CAS  PubMed  Google Scholar 

  38. Bilu C, Kronfeld-Schor N, Zimmet P, Einat H. Sex differences in the response to circadian disruption in diurnal sand rats. Chronobiol Int. 2022;39(2):169–85.

    CAS  PubMed  Google Scholar 

  39. Qian J, Morris CJ, Caputo R, Wang W, Garaulet M, Scheer F. Sex differences in the circadian misalignment effects on energy regulation. Proc Natl Acad Sci U S A. 2019;116(47):23806–12.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. Elfwing M, Nätt D, Goerlich-Jansson VC, Persson M, Hjelm J, Jensen P. Early stress causes sex-specific, life-long changes in behaviour, levels of gonadal hormones, and gene expression in chickens. PLoS One. 2015;10(5):e0125808.

    PubMed Central  Google Scholar 

  41. Bellizzi D, D’Aquila P, Scafone T, Giordano M, Riso V, Riccio A, et al. The control region of mitochondrial DNA shows an unusual CpG and non-CpG methylation pattern. DNA Res. 2013;20(6):537–47.

    CAS  PubMed Central  Google Scholar 

  42. Nedoluzhko A, Mjelle R, Renström M, Skjærven KH, Piferrer F, Fernandes JMO. The first mitochondrial 5-methylcytosine map in a non-model teleost (Oreochromis niloticus) reveals extensive strand-specific and non-CpG methylation. Genomics. 2021;113(5):3050–7.

    CAS  PubMed  Google Scholar 

  43. Patil V, Cuenin C, Chung F, Aguilera JRR, Fernandez-Jimenez N, Romero-Garmendia I, et al. Human mitochondrial DNA is extensively methylated in a non-CpG context. Nucleic Acids Res. 2019;47(19):10072–85.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. Lee JE, Oney M, Frizzell K, Phadnis N, Hollien J. Drosophila melanogaster activating transcription factor 4 regulates glycolysis during endoplasmic reticulum stress. G3 (Bethesda). 2015;5(4):667–75.

    PubMed  Google Scholar 

  45. Cagan A, Baez-Ortega A, Brzozowska N, Abascal F, Coorens THH, Sanders MA, et al. Somatic mutation rates scale with lifespan across mammals. Nature. 2022;604(7906):517–24.

  46. Sun JH, Ai SM, Liu SQ. Methylation-driven model for analysis of dinucleotide evolution in genomes. Theor Biol Med Model. 2020;17(1):3.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Kotrys AV, Szczesny RJ. Mitochondrial gene expression and beyond—novel aspects of cellular physiology. Cells. 2019;9(1):1–25.

  48. Black M, Arumugam P, Shukla S, Pradhan A, Ustiyan V, Milewski D, et al. FOXM1 nuclear transcription factor translocates into mitochondria and inhibits oxidative phosphorylation. Mol Biol Cell. 2020;31(13):1411–24.

    CAS  PubMed  Google Scholar 

  49. Macias E, Rao D, Carbajal S, Kiguchi K, Digiovanni J. Stat3 binds to mtDNA and regulates mitochondrial gene expression in keratinocytes. J Investig Dermatol. 2014;134(7):1971–80.

    CAS  PubMed  Google Scholar 

  50. Bird A. DNA methylation patterns and epigenetic memory. Genes Dev. 2002;16(1):6–21.

    CAS  PubMed  Google Scholar 

  51. Héberlé É, Bardet AF. Sensitivity of transcription factors to DNA methylation. Essays Biochem. 2019;63(6):727–41.

    PubMed  PubMed Central  Google Scholar 

  52. Barry WE, Thummel CS. The Drosophila HNF4 nuclear receptor promotes glucose-stimulated insulin secretion and mitochondrial function in adults. eLife. 2016;5:e11183.

    PubMed  PubMed Central  Google Scholar 

  53. Sasaki K, Uchiumi T, Toshima T, Yagi M, Do Y, Hirai H, et al. Mitochondrial translation inhibition triggers ATF4 activation, leading to integrated stress response but not to mitochondrial unfolded protein response. Biosci Rep. 2020;40(11):1–12.

    Google Scholar 

  54. Nargund AM, Fiorese CJ, Pellegrino MW, Deng P, Haynes CM. Mitochondrial and nuclear accumulation of the transcription factor ATFS-1 promotes OXPHOS recovery during the UPR(mt). Mol Cell. 2015;58(1):123–33.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. Quirós PM, Prado MA, Zamboni N, D’Amico D, Williams RW, Finley D, et al. Multi-omics analysis identifies ATF4 as a key regulator of the mitochondrial stress response in mammals. J Cell Biol. 2017;216(7):2027–45.

    PubMed  PubMed Central  Google Scholar 

  56. Landick R. Transcriptional Pausing as a Mediator of Bacterial Gene Regulation. Annu Rev Microbiol. 2021;75:291–314.

    CAS  PubMed  Google Scholar 

  57. Barshad G, Marom S, Cohen T, Mishmar D. Mitochondrial DNA Transcription and Its Regulation: An Evolutionary Perspective. Trends Genet. 2018;34(9):682–92.

    CAS  PubMed  Google Scholar 

  58. Yue Y, Ren L, Zhang C, Miao K, Tan K, Yang Q. Mitochondrial genome undergoes de novo DNA methylation that protects mtDNA against oxidative damage during the peri-implantation window. Proc Natl Acad Sci. 2022;119(30):1–9.

  59. Nitta KR, Jolma A, Yin Y, Morgunova E, Kivioja T, Akhtar J, et al. Conservation of transcription factor binding specificities across 600 million years of bilateria evolution. eLife. 2015;4:e04837.

    PubMed  Google Scholar 

  60. Palanker L, Tennessen JM, Lam G, Thummel CS. Drosophila HNF4 regulates lipid mobilization and beta-oxidation. Cell Metab. 2009;9(3):228–39.

    CAS  PubMed  PubMed Central  Google Scholar 

  61. Yang L, Chen Z, Stout ES, Delerue F, Ittner LM, Wilkins MR, et al. Methylation of a CGATA element inhibits binding and regulation by GATA-1. Nat Commun. 2020;11(1):2560.

    CAS  PubMed  PubMed Central  Google Scholar 

  62. Yang J, Horton JR, Wang D, Ren R, Li J, Sun D, et al. Structural basis for effects of CpA modifications on C/EBPβ binding of DNA. Nucleic Acids Res. 2019;47(4):1774–85.

    CAS  PubMed  Google Scholar 

  63. Hashimoto H, Wang D, Horton JR, Zhang X, Corces VG, Cheng X. Structural Basis for the Versatile and Methylation-Dependent Binding of CTCF to DNA. Mol Cell. 2017;66(5):711-20.e3.

    CAS  PubMed  PubMed Central  Google Scholar 

  64. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009;462(7271):315–22.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. Ziller MJ, Müller F, Liao J, Zhang Y, Gu H, Bock C, et al. Genomic distribution and inter-sample variation of non-CpG methylation across human cell types. PLoS Genet. 2011;7(12):e1002389.

    CAS  PubMed  PubMed Central  Google Scholar 

  66. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    PubMed  PubMed Central  Google Scholar 

  67. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    CAS  PubMed  Google Scholar 

  68. Robinson MD, Stirzaker C, Statham AL, Coolen MW, Song JZ, Nair SS, et al. Evaluation of affinity-based genome-wide DNA methylation data: effects of CpG density, amplification bias, and copy number variation. Genome Res. 2010;20(12):1719–29.

    CAS  PubMed  Google Scholar 

  69. Grant CE, Bailey TL. XSTREME: Comprehensive motif analysis of biological sequence datasets. bioRxiv. 2021;458722:1–3.

Download references

Acknowledgements

The authors have no additional acknowledgements.

Funding

Open access funding provided by Uppsala University. The authors would like to thank the following funding agencies for their contributions to the research presented:

Formas 2019–02084 (PJ) and 2018–01074 (CGB).

European Research Council Advanced grant 322,206 (“Genewell”) (PJ).

The Swedish Research council (VR) 2019-04053_VR (CGB, JL).

The São Paulo Research Foundation (FAPESP) projects #2016/20440–3 and #2018/13600–0 (FP).

Author information

Author notes

  1. John Lees and Fábio Pèrtille contributed equally to this work.

Authors and Affiliations

  1. Evolutionsbiologiskt Centrum (EBC), Uppsala University, Uppsala, 75236, Sweden

    John Lees, Fábio Pèrtille & Carlos Guerrero Bosagna

  2. Institutionen För Fysik, Kemi Och Biologi (IFM), Linköping University, Linköping, 58330, Sweden

    Pia Løtvedt & Per Jensen

Authors
  1. John Lees
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Fábio Pèrtille
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Pia Løtvedt
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Per Jensen
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Carlos Guerrero Bosagna
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Experimental design: PJ. Experimental work and sample collection: PL. Conceptualization: JL, FP, CGB. Investigation: JL, FP, CGB. Visualization: JL, FP. Supervision: CGB, PJ. Writing—original draft: JL, FP. Writing—review & editing: PJ, CGP, FP, PL, JL.

Corresponding author

Correspondence to Carlos Guerrero Bosagna.

Ethics declarations

Ethics approval and consent to participate

Experimental protocols were approved by the Regional Council for Ethical Licensing of Animal Experiments (Linköpings djurförsöksetiska nämnd) and conducted in accordance with all relevant guidelines and regulations (License Nr. 50–13 Linköping, Sweden). Where relevant to the study, all methods are reported in accordance with the ARRIVE guidelines.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1.

Locations of DNA-binding motifs for HNF4 (a) and ZNF324 (b) (taken from a human database and locations identified using FIMO) within the mitochondrial genomes of Gallus gallus, Homo sapiens, Alligator mississippiensis, Danio rerio and Drosophila melanogaster. Fig. S2. Schematic of the bioinformatics pipeline used to identify the differentially methylated based in a CpN context throughout the mitochondrial genome. Table S1. Table of the raw sequencing read data, corresponding to the whole genome (WG), the mitochondrial chromosome (MT) and representing true mitochondrial genome reads after the removal of nuclear mitochondrial pseudogenes (numts). Table S2. Start base, stop base, length, gene location within the mtDNA, mean direction of cytosine methylation and sequence of the mtDMRs identified in our analysis across different pairwise statistical comparisons. Contrasts representing statistical comparisons are as follows; male and female control birds (MC vs. FC), male and female stressed birds (MS vs. FS), male control vs. male stress birds (MC vs. MS) and female control vs. female stress birds (FC vs. FS). Table S3. Base composition within mtDMRs across our pairwise statistical comparisons compared to the base composition of the mitochondrial genome (A = 30.3%, T = 23.8%, C = 32.5%, G = 13.5%). Contrasts representing statistical comparisons are as follows; male and female control birds (MC vs. FC), male and female stressed birds (MS vs. FS), male control vs. male stress birds (MC vs. MS) and female control vs. female stress birds (FC vs. FS). Table S4. A comparison of the occurrences CpN dinucleotides within mtDMRs. CpA, CpT, CpC and CpC occurrences within mtDMRs of the mitcochondrial DNA of male and female chickens exposed to circadian light stress. Values are expressed in relation to pairwise statistical comparisons between Male control vs. female control (MC vs. FC), male stress vs. female stress (MS vs. FS), male control vs. male stress (MC vs. MS) and female control vs. female stress (FC vs. FS). Chi squared p-values > 0.05 indicate no difference between the number of CpNs with differential sequencing coverage in comparison to the total number of CpNs present within the mtDMRs. Table S5. Frequency matrix for human ATF4 used to determine ATF4 locations within diverse taxa. Source:  http://floresta.eead.csic.es/footprintdb/index.php?motif=591111e9f6326e96b78d40f0820b6e8d. Table S6. Frequency matrix for human ZNF324 used to determine ZNF324 locations within diverse taxa. Source: https://jaspar2022.genereg.net/matrix/MA1977.1/. Table S7. Frequency matrix for human HNF4A used to determine HNF4A locations within diverse taxa. Source: https://jaspar2022.genereg.net/matrix/MA0114.4/. Table S8. Locations of DNA-binding motifs for ATF4 (taken from a human database and locations identified using FIMO) within the mitochondrial genomes of Gallus gallus, Homo sapiens, Alligator mississippiensis, Danio rerio and Drosophila melanogaster. Table S9. Locations of DNA-binding motifs for ZNF324 (taken from a human database and locations identified using FIMO) within the mitochondrial genomes of Gallus gallus, Homo sapiens, Alligator mississippiensis, Danio rerio and Drosophila melanogaster. Table S10. Locations of DNA-binding motifs for HNF4 (taken from a human database and locations identified using FIMO) within the mitochondrial genomes of Gallus gallus, Homo sapiens, Alligator mississippiensis, Danio rerio and Drosophila melanogaster. Table S11. Table used to calculate the CpG enrichment score of samples based upon MeDIP samples in comparison to the reference genome. Data S1. The differentially methylated CpN dinucleotides identified through pairwise statistical comparisons and their locations within the mitochondrial genome. Contrasts representing statistical comparisons are as follows; male and female control birds (MC vs. FC), male and female stressed birds (MS vs. FS), male control vs. male stress birds (MC vs. MS) and female control vs. female stress birds (FC vs. FS).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lees, J., Pèrtille, F., Løtvedt, P. et al. The mitoepigenome responds to stress, suggesting novel mito-nuclear interactions in vertebrates. BMC Genomics 24, 561 (2023). https://doi.org/10.1186/s12864-023-09668-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09668-9

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us