High-resolution profiles of gene expression and DNA methylation ...

Advance Publication by J-STAGE Journal of Reproduction and Development

Accepted for publication: February 16, 2017 Published online: March 31, 2017

1

High-resolution profiles of gene expression and DNA methylation highlight mitochondrial

2

modifications during early embryonic development

3

Likun Ren1)*, Chao Zhang1)*, Li Tao1)*, Jing Hao1), Kun Tan1), Kai Miao1), Yong Yu1), Linlin Sui1),

4

Zhonghong Wu1), Jianhui Tian 1), Lei An1)

5

1) Key Laboratory of Animal Genetics, Breeding and Reproduction of the Ministry of Agriculture,

6

National Engineering Laboratory for Animal Breeding, College of Animal Science and Technology,

7

China Agricultural University, Beijing 100193, P. R. China

8

* Likun Ren and Chao Zhang contributed equally to this work.

9

Running head: Mitochondrial modifications in embryos

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

Corresponding author: Lei An (email: [email protected])

26

Abstract

27

Well-organized mitochondrial functions and dynamics are critical for early embryonic development

28

and are operated via a large number of mitochondria-related genes (MtGs) encoded by both the

29

nuclear and the mitochondrial genome. However, the mechanisms underlying mitochondrial

30

modifications during the critical window between blastocyst implantation and postimplantation

31

organogenesis are poorly understood. Herein, we performed high-resolution dynamic profiling of

32

MtGs to acquire a more detailed understanding of mitochondrial modifications during early

33

development. Our data suggest that the resumption of mitochondrial mass growth is not only

34

facilitated by increased mitochondrial biogenesis and mitochondrial DNA (mtDNA) replication, but

35

also by the appropriate balance between mitochondrial fission and fusion. In addition, increased levels

36

of reactive oxygen species (ROS) resulting from enhanced mitochondrial functions may be the critical

37

inducer for activating the glutathione (GSH)-based stress response system in early embryos. The

38

appropriate balance between the mitochondrial stress response and apoptosis appears to be significant

39

for cell differentiation and early organogenesis. Furthermore, we found that most MtGs undergo de

40

novo promoter methylation, which may have functional consequences on mitochondrial functions and

41

dynamics during early development. We also report that mtDNA methylation can be observed as early

42

as soon after implantation. DNMT1, the predominant enzyme for maintaining DNA methylation,

43

localized to the mitochondria and bound to mtDNA by the implantation stage. Our study provides a

44

new insight into the involvement of mitochondria in early mammalian embryogenesis. We also

45

propose that the epigenetic modifications during early development are significant for modulating

46

mitochondrial functions and dynamics.

47

Key words: Mitochondria, Early embryos, Reactive oxygen species, Glutathione, DNA methylation

48

49

Introduction

50

Mitochondria not only provide the energy required to maintain cellular survival and growth but

51

also participate in a number of pathways that maintain cellular homeostasis, including ion homeostasis,

52

amino acid metabolism, glycolysis, fatty acid metabolism, signal transduction, and apoptosis [1].

53

Well-organized mitochondrial dynamics and functions are of prime importance for early embryonic

54

development. During the period between the blastocyst and the postimplantation stages, which

55

overlaps with implantation and organogenesis initiation, the mitochondrion undergoes dramatic

56

modifications, including mitochondrial biogenesis and mitochondrial DNA (mtDNA) replication, a

57

shift from anaerobic to aerobic metabolism, and changes in the mitochondrial ultrastructure [2]. These

58

mitochondrial modifications are essential for successful implantation and survival. However, the

59

comprehensive mechanism underlying mitochondrial dynamics during this critical developmental

60

period remains largely elusive.

61

The orchestrated mitochondrial functions and dynamics are highly regulated by both the nuclear

62

(nDNA) and mitochondrial genomes (nuclear-mitochondrial cross-talk). In addition to the 37

63

mtDNA-encoded genes, it is estimated that approximately 1200 proteins encoded by nDNA can enter

64

the mouse mitochondrial matrix or bind to its membrane [3]. These mitochondria-related genes (MtGs)

65

have been targeted in high-resolution gene expression analysis to evaluate the involvement of

66

mitochondrial dysfunctions in cardiac pathology [4] and cancer [5]. Recently, our own studies

67

suggested that mitochondrial dysfunction may be the determining factor of impaired development of

68

in vitro fertilized (IVF) embryos through disturbed oxidative phosphorylation (OXPHOS), reduced

69

mitochondrial biogenesis, and dysregulated responses to oxidative stress [6].

70

In this study, we performed a dynamic high-resolution expression profiling of MtGs during early

71

embryonic development, using samples from E3.5 blastocysts, epiblasts from E7.5 embryos (E7.5

72

epiblasts), and E10.5 embryos. The targeted profiling of a specific gene set leads to a more efficient

73

enrichment of the involved mitochondrial functions, thereby providing a more detailed and

74

comprehensive understanding of mitochondrial dynamics and functions during early development.

75

In addition to the expression profiling, we also highlighted the DNA methylation dynamics of

76

MtGs for the period between the blastocyst and the postimplantation stages. One of the most important

77

epigenetic events during this critical developmental period is the re-establishment of DNA methylation

78

patterns, termed de novo DNA methylation. The progression of this. epigenetic event highly coincides

79

with the mitochondrial modifications taking place by the implantation stage, such as the ‘embryonic

80

shift’ from anaerobic to aerobic metabolism and the resumption of mitochondrial biogenesis. In

81

addition, previous studies demonstrated that the DNA methylation status of nDNA-encoded MtGs

82

underlies tissue- or cell type-dependent mitochondrial functions and dynamics [7, 8] as well as

83

mitochondrial pathology [9-11]. These facts lead us to investigate whether de novo DNA methylation

84

may play a regulatory role in the mitochondria during early development.

85

Materials and Methods

86

Animal preparation

87

F1 female mice (ICR, 5–6 weeks old) and F1 male mice (ICR, 8–9 weeks old) were fed ad

88

libitum and housed in a room with a controlled light cycle (12L:12D). F1 females were naturally

89

mated with F1 males. All studies were specifically approved by and performed in accordance with the

90

guidelines of the Institutional Animal Care and Use Committee of China Agricultural University.

91

Embryo preparation and collection

92

The criteria for sampling embryos were based on developmental progress and morphology.

93

Embryos showing typical morphological features according to the well-established landmarks [12]

94

were sampled for further analyses. All sampled embryos for RNA and DNA isolation were serially

95

washed with phosphate-buffered saline (PBS; GIBCO, Life Technologies, NY, USA) and immediately

96

stored in liquid nitrogen for future use. In addition, to avoid cross-contamination between embryonic

97

and extraembryonic tissues, we controlled the efficiency of the dissection procedure by detecting the

98

expression of markers specific to the extraembryonic ectoderm (ETS-related transcription factor, Elf5),

99

ectoplacental cone (Achaete-scute like 2, ASCL2, also known as Mash2), and epiblasts (fibroblast

100

growth factor 5, Fgf5) at E7.5 [13].

101

In vitro fertilization and culture

102

To generate in vitro fertilized embryos, female ICR mice were superovulated by an

103

intraperitoneal (i.p.) injection of 5 IU pregnant mare serum gonadotropin (PMSG; Ningbo Hormone

104

Product Co., Ltd., Ningbo, China), followed by an i.p. injection of 5 IU human chorionic gonadotropin

105

(hCG; Ningbo Hormone Product Co., Ltd.) 48 h later. At 14 h after hCG treatment, cumulus-enclosed

106

oocyte complexes were recovered from oviducts and cumulus cells were removed by digestion with

107

hyaluronidase (Sigma-Aldrich, St. Louis, MO, USA) for 3–5 min. The oocytes were rinsed in human

108

tubal fluid (HTF) medium (Sage, Bedminster, NJ, USA), and placed into 60 μl drops of HTF medium

109

covered with paraffin oil, before being equilibrated overnight at 37°C and 5% CO2. Sperm was

110

collected from the cauda epididymis and capacitated for 1 h in HTF medium at 37°C and 5% CO2.

111

Oocytes were inseminated 15 h post-hCG with 106 spermatozoa. After 4 h in the incubator, oocytes

112

and zygotes were washed several times in potassium simplex optimization medium containing amino

113

acids (KSOM + AA; Millipore, Billerica, MA, USA) and then transferred to 60 µl drops of KSOM

114

medium. The zygotes, determined by the presence of two pronuclei, were cultured to the blastocyst

115

stage at 37°C in a 5% CO2 atmosphere. Well-developed late-cavitating blastocysts of similar

116

morphology were sampled at 106–112 h post-hCG after culturing in KSOM medium [14].

117

Pseudopregnant female mice were mated with vasectomized males 3.5 days prior to embryo transfer.

118

Twelve IVF blastocysts were transferred to a single recipient, with six embryos placed in each uterine

119

horn.

120

High-throughput RNA sequencing (RNA-seq) and biological analysis

121

Total RNAs were extracted from embryos at different stages with the TRIzol Reagent (Invitrogen,

122

Carlsbad, CA, USA). Polyadenylated RNAs were isolated using the Oligotex mRNA Midi Kit (Qiagen,

123

Valencia, CA, USA). The RNA-seq libraries were constructed using the SOLiD Whole Transcriptome

124

Analysis Kit (Applied Biosystems, Carlsbad, CA, USA) following the standard protocol and

125

sequenced on the SOLiD platform (Applied Biosystems) to generate high-quality single-end reads.

126

The raw reads were aligned to genome sequences, trimming off a nucleotide from each of the 5′ and 3′

127

ends and allowing up to two mismatches. Reads mapped to multiple locations were discarded and only

128

uniquely mapped reads were used for the subsequent analysis. Gene expression levels were measured

129

in reads per kilobase of exon model per million mapped reads (RPKM). A single pooling strategy for

130

relatively large sample sizes was used, as described previously [15, 16].

131

DAVID v6.7 (http://david.abcc.ncifcrf.gov) was used to annotate biological themes (gene

132

ontology, GO). REVIGO analysis (http://revigo.irb.hr) was also performed to visualize the interactive

133

relationship among enriched processes. The KEGG database (http://www.genome.jp/kegg) was used to

134

determine the associated pathways. Data on MtGs were fed into STRING (http://string.embl.de) to

135

build an interaction network. The phenotype annotations of MtGs were analyzed based on the Mouse

136

Genome Informatics (MGI; http://www.informatics.jax.org/phenotypes.shtml) database.

137

Apoptosis analysis

138

Collected blastocysts were washed three times with 0.1% PVA/PBS and then transferred to PBS

139

supplemented with 4% (v/v) paraformaldehyde and 0.5% Triton X-100 for simultaneous fixation and

140

permeabilization at 37°C for 45 min. Apoptotic nuclei were detected using the In Situ Cell Death

141

Detection Kit (Roche, Mannheim, Germany). The number of apoptotic nuclei was determined under a

142

BX51epifluorescence microscope (Olympus, Tokyo, Japan).

143

RNA isolation and quantitative real-time RT-PCR (qRT-PCR)

144

Total RNA extraction from embryos was performed using TRIzol according to the manufacturer’s

145

instructions. Reverse transcription (RT) was carried out using a commercially available first-strand

146

cDNA synthesis kit (iScript cDNA Synthesis Kit; Bio-Rad Laboratories, USA). Real-time PCR was

147

performed in the Bio-Rad CFX96 Real-Time PCR System using SsoFast EvaGreen Supermix

148

(Bio-Rad Laboratories, Hercules, CA, USA). Used primers are listed in Supplementary Table S1.

149

MtDNA copy number determination by qPCR

150

The nuclear-encoded β-actin gene and the mitochondria-encoded ND5 gene were used to

151

determine the mtDNA copy numbers with the primers listed in Supplementary Table S1. Standard

152

curves were constructed using recombinant plasmids that were serially diluted 103–108 times. The

153

mtDNA and nDNA copy numbers were calculated from the threshold cycle number (CT) and

154

corrected using the standard curve. The overall number of mtDNA copies per cell was calculated using

155

the following formula: mtDNA copies per cell = ND5 copy number / (β-actin copy number/2).

156

Measurement of reduced glutathione (GSH) and oxidized glutathione (GSSG)

157

The contents of intracellular reduced glutathione (GSH) and oxidized glutathione (GSSG) in each

158

blastocyst were determined using commercially available kits (Beyotime, Jiangsu, China) according to

159

the manufacturer’s instructions. Absorbance was measured at 412 nm using the Infinite M200 PRO

160

NanoQuant microplate reader (Tecan, Männedorf, Switzerland). Approximately 30 E3.5 blastocysts

161

and 20 E7.5 epiblasts were used for measurements.

162

Methylated DNA immunoprecipitation-sequencing (MeDIP-seq)

163

Global genomic DNA, including both nDNA and mtDNA, was isolated from embryonic samples

164

using the DNeasy Kit (Qiagen) according to the manufacturer’s instructions. The quality of each DNA

165

sample with respect to integrity, purity, and concentration was assessed using a DN-1000

166

Spectrophotometer (Nanodrop Technologies, Wilmington, DE, USA). Genomic DNA was then

167

fragmented using a Covarias sonication system (Covarias, Woburn, MA, USA). After sonication, the

168

fragments were denatured to produce single-stranded DNA (ssDNA). Following denaturation, the

169

ssDNA was incubated with antibodies recognizing 5-methylcytosine (5mC). Magnetic beads

170

conjugated to anti-mouse-IgG were then used to bind the anti-5mC antibodies and unbound DNA was

171

removed along with the supernatant. Finally, DNA was released by digesting the antibodies with

172

proteinase K and collected. Sequencing was carried out in 50 bp reads on the Illumina Hiseq 2000

173

(Illumina, San Diego, CA, USA) by the Beijing Genomics Institute (BGI, Shenzhen, Guangdong,

174

China) following the standard protocol.

175

MeDIP-qPCR

176

Global DNA was isolated using the Solar MicroElute Genomic DNA Kit (Solar Technologies Inc.,

177

Gaithersburg, MD, USA). Sonicated DNA (fragment length: 100–600 bp) was treated with the

178

EpiQuik MeDIP Ultra Kit (Epigentek, Farmingdale, NY, USA) to enrich methylated DNA fragments.

179

Briefly, samples were precipitated with anti-5-mC, while samples precipitated with non-immunized

180

IgG served as a negative control. The enriched methylated DNA was quantified by qPCR using

181

specific primers (Supplementary Table S1) and the enrichment efficiency of MeDIP-qPCR was

182

measured as fold enrichment (FE) with the following formula: FE (%) = 2(IgG CT –Sample CT) × 100%, as

183

presented in the product manual.

184

Results

185

General dynamic profiling of MtGs during early development

186

We gathered information for over 1000 MtGs, as reported by Mootha’s study [3], from the

187

transcriptomic data obtained from E3.5 blastocysts, E7.5 epiblasts, and E10.5 embryos

188

(Supplementary Table S2). The hierarchical clustering heatmap showed that all MtGs can be sorted

189

into three main groups based on expression dynamics, namely whether they were most abundant at

190

E3.5, E7.5, or E10.5. Each group can be further divided into two clusters (Fig. 1A). Compared with

191

E3.5 blastocysts, E7.5 epiblasts and E10.5 embryos were clustered more tightly, with similar MtG

192

expression patterns. Functional analysis showed that each cluster was significantly associated with

193

many aspects of mitochondrial functions or related processes, including OXPHOS, oxoacid and amino

194

acid metabolism, and apoptosis (Supplementary Fig. 1A), implying that mitochondrial dynamics

195

during early development are orchestrated via complicated expression patterns involving MtGs

196

(Supplementary Fig. 1B).

197

To further explore the involvement of mitochondrial functions and dynamics during early

198

embryonic development, the expression patterns of MtGs at E3.5, E7.5, and E10.5 were compared

199

consecutively, i.e., E7.5 versus E3.5 and E10.5 versus E7.5. These two transitions encompass a series

200

of important developmental events, such as implantation, cellular differentiation, and organogenesis.

201

In the first transition, i.e., from E3.5 blastocysts to E7.5 epiblasts, 920 MtGs displayed significant

202

changes in expression (P < 0.05), accounting for approximately 80% of identified MtGs. During the

203

transition from E7.5 epiblasts to E10.5 embryos, the expression of 843 MtGs changed significantly (P

204

< 0.05), accounting for approximately 73% of identified MtGs. Both percentages were higher than the

205

corresponding values of the global transcriptome, which were approximately 70% and 64%,

206

respectively (Fig. 1B), implying that expression of MtGs is highly dynamic during early development.

207

These observations also suggest that mitochondrial features may undergo a more dramatic transition

208

from E3.5 blastocysts to E7.5 epiblasts, in accordance with the results of the hierarchical clustering.

209

A Venn diagram analysis was performed to screen MtGs consistently up- or downregulated (FC >

210

2) during the transition E3.5 to E10.5. According to the results of the analysis, 58 MtGs showed

211

consistent upregulation while 31 MtGs were consistently downregulated (Fig. 1C). Considering that

212

the expression of these MtGs consistently changed during early development, they may contribute

213

significantly to the highly dynamic mitochondrial functions taking place during this stage.

214

Detailed functional profiling of MtGs during early development

215

To further understand the comprehensive mechanism underlying mitochondrial functions and

216

dynamics during early development, high-resolution functional profiling was performed for the MtGs

217

whose expression significantly changed (P < 0.05, FC > 2) during the transition from E3.5 to E7.5

218

(Supplementary Table S3) as well as from 7.5 to E10.5 (Supplementary Table S4).

219

GO classification of biological processes indicated that ‘Oxoacid metabolic process’, ‘Coenzyme

220

metabolic process’, ‘Cellular amine metabolic process’, and ‘Fatty acid metabolic process’ were the

221

four most significantly enriched terms in the transition from E3.5 blastocysts to E7.5 epiblasts,

222

suggesting that metabolic changes are the most evident feature of mitochondrial modifications during

223

this period. In addition, ‘Apoptotic mitochondrial changes’ and ‘Response to ROS’ may also be

224

significantly involved in this transition. Similar processes were enriched in the transition from E7.5

225

epiblasts to E10.5 embryos (Fig. 2A and B). Next, we performed REVIGO analysis to visualize the

226

interactive relationship among the enriched processes. Results from both transitions suggest that

227

‘Transmembrane transport’ may play a key role in regulating mitochondrial metabolism and

228

mitochondria-mediated apoptosis (Fig. 2C and D).

229

KEGG pathway analysis provided us with more detailed insights into mitochondrial

230

modifications. For the transition from E3.5 blastocysts to E7.5 epiblasts, the altered mitochondrial

231

energy metabolism included ‘TCA cycle’, ‘OXPHOS’, and ‘Fatty acid β-oxidation’, while the changed

232

amino acid metabolism involved arginine, proline, valine, leucine, glutamate, and glycine. It should be

233

noted that some pathways of neurodegenerative diseases, including ‘Alzheimer’s disease (AD)’,

234

‘Parkinson’s disease (PD)’ and ‘Huntington’s disease (HD)’ were also enriched, mainly because of

235

MtGs associated with ‘OXPHOS’ and ‘apoptosis’. Likewise, changes in energy metabolism and amino

236

acid metabolism were also significant during the transition from E7.5 to E10.5 embryos (Table 1).

237

Using dynamically changed MtGs (P < 0.05, FC > 2) as seed nodes, interaction networks were

238

constructed. In addition to functional clusters similar to those identified by the GO and KEGG

239

analyses, other processes or functions, such as ‘mitochondrial translation’ and ‘mitochondrial

240

biogenesis’, were also significantly clustered. Another interesting observation was that some

241

well-known genes responsible for the regulation of cell pluripotency and differentiation, such as

242

members of the MAPK (Mapk8 and Mapk12), WNT (Wnt1 and Wnt3), and SMAD (Smad3 and Smad7)

243

families, were clustered tightly with other MtGs, suggesting an essential role of mitochondrial

244

functions in the regulation of pluripotency and differentiation in early embryos (Fig. 2E and F).

245

Mitochondrial modifications required for successful implantation

246

The above-mentioned functional profiling revealed that many mitochondrial modifications had

247

undergone evident temporal changes by the implantation stage. The observed changes in the

248

expression of some key genes of energy metabolism were in accordance with the ‘embryonic shift’

249

from anaerobic to aerobic metabolism. Accordingly, MtGs associated with mitochondrial biogenesis

250

showed changes that may accommodate the enhanced aerobic metabolism. It is noteworthy that the

251

expression of MtGs underlying GSH biosynthesis and functions changed correspondingly, suggesting

252

a physiological response to the increased ROS levels resulting from the enhanced mitochondrial

253

functions (Fig. 3A).

254

To gain further insight into the comprehensive mechanism underlying the increased

255

mitochondrial biogenesis and GSH production in early embryonic development, we analyzed the

256

dynamics of MtGs involved in these processes based on the gene list provided by Quick GO

257

(https://www.ebi.ac.uk/QuickGO)

and

Wikipathways

258

(http://www.wikipathways.org/index.php/WikiPathways). Results showed that the expression of most

259

of the MtGs associated with these processes correspondingly changed from E3.5 to E7.5 (Fig. 3B).

260

The changes of some selected MtGs were validated using qRT-PCR (Fig. 3C). Next, to understand the

261

interactive relationship of these processes, we constructed a small network. This revealed a high level

262

of association among these processes (Fig. 3D). Indeed, the involvement of mitochondrial biogenesis

263

and GSH-based stress response in mitochondria-mediated apoptosis can also be deduced by the fact

264

that many MtGs are shared among these terms.

265

To assess the association between the resumption of mitochondrial biogenesis and mtDNA

266

replication on the embryonic developmental potential, we compared the mtDNA copy numbers

267

between in vivo conceived (IVO) and in vitro fertilized (IVF) embryos, which are characterized by

268

higher and lower developmental rates, respectively. Our previous studies showed that IVF embryos

269

have lower survival rate and weight, as well as greater variability in size and increased growth defects,

270

compared to their IVO counterparts [6, 13, 17]. Quantitative analysis showed that mtDNA copy

271

numbers in IVF embryos laid in a significantly lower range compared to those in IVO embryos (Fig.

272

3E), suggesting that depressed mitochondrial resumption is strongly associated with impaired

273

embryonic developmental potential.

274

The presence of enhanced mitochondrial mass and metabolism during the implantation stage led

275

us to speculate that postimplantation embryos face increased oxidative stress compared to

276

preimplantation embryos. To this end, we determined the GSH/GSSG ratio, which is a well-accepted

277

indicator for evaluating oxidative stress status [18], in E3.5 blastocysts and E7.5 epiblasts. We found

278

the latter to display a significantly higher GSH/GSSG ratio, indicating increased oxidative stress soon

279

after implantation (Fig. 3F).

280

Subsequently, because of the significant changes in the expression patterns of apoptotic MtGs

281

that had taken place by the implantation stage, we decided to compare the occurrence of apoptosis in

282

E3.5 blastocysts and E7.5 epiblasts. TUNEL analysis showed that E7.5 epiblasts displayed a high rate

283

of apoptosis compared to E3.5 blastocysts (Fig. 3G). Possible false-negative apoptotic nuclei in E3.5

284

blastocysts can be excluded by the parallel observation of TUNEL-positive apoptotic nuclei in IVF

285

blastocysts, which have been reported to have a relative higher apoptotic rate [6].

286

The involvement of MtGs in embryonic survival and organogenesis

287

To test the potential role of MtGs in embryonic survival and organogenesis, we performed an

288

MGI analysis using dynamically changed MtGs (P < 0.05, FC >2; Supplementary Tables S5 and S6).

289

Many MtGs were annotated with phenotypes associated with aberrant embryonic/fetal development,

290

such as ‘Complete embryonic lethality’, ‘Embryonic growth arrest’, ‘Failure of somite differentiation’,

291

and ‘Failure of embryo implantation’. In addition, many MtGs were associated with abnormal

292

morphology or aberrant development of organs, especially of the neurosensory and cardiovascular

293

systems, as well as with lethality during organogenesis. These observations suggested that orchestrated

294

mitochondrial functions and dynamics are essential for embryonic survival and normal organogenesis

295

both during and after implantation. In addition, we found that different mitochondrial modifications

296

might functionally correlate with each other to regulate early organogenesis (Fig. 4).

297

Dynamics of promoter DNA methylation of nDNA-encoded MtGs during early development

298

The dramatic changes of mitochondrial functions and dynamics during early development,

299

including the resumption of mitochondrial biogenesis and the shift from anaerobic to aerobic

300

metabolism, coincide with de novo DNA methylation, which is one of the most important epigenetic

301

events during early development. This led us to hypothesize that the de novo DNA methylation of

302

MtGs may have an essential regulatory role in the mitochondrial modifications taking place during

303

early development. To test this hypothesis, nDNA-encoded MtGs whose promoters were methylated at

304

least once at three time-points (Supplementary Table S7) were hierarchically clustered based on their

305

promoter DNA methylation levels, as described previously [19]. The resulting heatmap showed a

306

generally progressive increase in promoter DNA methylation of MtGs during early development,

307

especially of those belonging to clusters 1 and 2, which correspond to genes that are unmethylated in

308

E3.5 blastocysts and gain promoter DNA methylation before E7.5 or E10.5, respectively. On the other

309

hand, promoter DNA methylation increased slightly in cluster 3, which includes genes whose

310

promoters appear methylated throughout early development (Fig. 5A).

311

To test the correlation between promoter de novo methylation and the expression levels of MtGs

312

during early development, we integrated the promoter DNA methylation profiles to the expression

313

profiles of MtGs in clusters 1 and 2 (Fig. 5B; boxed by green and blue dotted lines, respectively). In

314

cluster 1, 48 MtGs (54.5%) were downregulated during the transition from E7.5 to E10.5, with

315

increased levels of promoter DNA methylation. In cluster 2, 85 MtGs (45.5%) showed decreased

316

expression at E7.5 compared to E3.5, and the reduction in expression coincided with promoter de novo

317

methylation taking place before E7.5. This result suggests that promoter de novo methylation may

318

have a considerable regulatory effect on the expression of approximately half of MtGs, during early

319

development.

320

To further analyze the functional consequences of de novo methylation on mitochondrial

321

modifications during early development, GO and KEGG analyses were performed for MtGs that were

322

presumably regulated by increased promoter DNA methylation (indicated by green dots in Fig. 5A).

323

The results showed that these MtGs were significantly associated with most mitochondrial functions

324

taking place during both transitions. Amongst others, these included ‘Oxoacid metabolic process’,

325

‘Cellular amino acid metabolic process’, ‘Electron transport chain’, ‘Mitochondrial transport’, and

326

‘OXPHOS’ (Fig. 5C and D). It should be noted that although the number of enriched MtGs was very

327

low (2–7), the levels of significance were very high, which implied that MtG promoter de novo

328

methylation may profoundly affect mitochondrial functions and dynamics during early development.

329

MtDNA methylation in postimplantation embryos

330

Dnmt1 expression was found to be significantly higher in postimplantation than preimplantation

331

embryos (Fig. 6A). This observation was confirmed by qRT-PCR analysis. As Dnmt1 can also

332

translocate to the mitochondria for maintaining mtDNA methylation [20], we further investigated the

333

expression dynamics of a transcript variant of Dnmt1 specifically targeted to mitochondria (mt-Dnmt1),

334

which is not fully synchronized with the levels of total Dnmt1 (Fig. 6B). Based on the relatively stable

335

expression of mt-Dnmt1 up to the implantation stage, we hypothesized that embryonic mtDNA may

336

undergo methylation after implantation. This idea was supported by the enrichment of MeDIP reads on

337

mtDNA in E7.5 epiblasts and E10.5 embryos (Fig. 6C, Supplementary Table S8). We validated this

338

observation by performing MeDIP-qPCR with a highly sensitive and specific anti-5mC antibody. It

339

should be mentioned that the primers used for detection were blasted against the mouse nuclear

340

genome to prevent a possible contamination from nuclear sequences of mitochondrial origin.

341

Focusing on the MeDIP-qPCR results from both hyper- (ND5, ND6) and hypomethylated (12S

342

rRNA) regions, as well as a control region (D-loop), we observed that each of these regions was

343

significantly methylated, with an increasing trend during the transition from E7.5 to E10.5 (Fig. 6D).

344

Moreover, to confirm that mt-Dnmt1 can translocate to the mitochondria by the implantation stage, we

345

tested the binding of mt-Dnmt1 to mtDNA by CHIP-qPCR using a non-specific anti-DNMT1 antibody

346

and primers specific to mtDNA in E6.5 embryos, which represent a stage immediately after

347

implantation. Results showed a significant enrichment of DNMT1 at two selected hypermethylated

348

regions (ND5 and ND6) of mtDNA (Fig. 6E).

349 350

Discussion

351

It is well documented that mitochondrial functions and dynamics are of prime importance to early

352

development. Well-organized mitochondrial modifications are thought to be essential for embryonic

353

survival and growth during implantation, as well as for subsequent differentiation and organogenesis.

354

Focusing on MtGs, this study provided a comprehensive and detailed understanding of the

355

mechanisms responsible for the orchestrated mitochondrial modifications taking place during early

356

development.

357

During the transition from E3.5 to E7.5, embryos undergo the differentiation of embryonic germ

358

layers and implantation, two of the most important challenges that enable further embryonic survival

359

and growth. On the other hand, in the period between E7.5 and E10.5, embryos initiate organogenesis

360

following the completion of gastrulation [21]. During the transition from the blastocyst to the

361

gastrulation stage, embryos undergo a shift in energy metabolism, in which glycolytic activity

362

gradually decreases but OXPHOS is rapidly enhanced. Coincident with this shift is the resumption of

363

mitochondrial biogenesis and mtDNA replication. Increased mitochondrial biogenesis and mtDNA

364

replication in somatic cells are thought to accommodate high energy demands during physiological

365

and pathological processes [22, 23]. In our study, we found that the expression of many key MtGs

366

involved in mitochondrial biogenesis and mtDNA replication significantly increases during

367

implantation. Among these, some are well-established MtGs considered responsible for the resumption

368

of mitochondrial biogenesis in early embryos, such as Polg, Polg2, and Tfam [24-26]. Other MtGs that

369

have been previously reported to play important roles in somatic cells may also be critical for this

370

resumption. For example, Ssbp1 encodes a subunit of a single-stranded DNA-binding complex that

371

can maintain the stability of mtDNA by removing mismatches during mtDNA replication [27] and

372

support mitochondrial protein synthesis and morphology in human osteosarcoma and embryonic

373

kidney cells [28, 29]. Stoml2, also known as Slp2, codes for a widely expressed mitochondrial inner

374

membrane protein and can induce mitochondrial biogenesis in human T lymphocytes [30]. Tk2 codes

375

for a deoxyribonucleoside kinase that localizes to the mitochondria and is required for the replication

376

and maintenance of mtDNA in cultured human somatic cells [31]. A recent study demonstrated that

377

Tk2 may be involved in the early development of zebrafish by modulating mtDNA metabolism [32];

378

however, a role for this gene in the resumption of mtDNA replication in mammalian embryos has not

379

yet been identified.

380

Mitochondria are dynamic organelles that are constantly undergoing fission and fusion to

381

maintain mitochondrial homeostasis and morphology [33-36]. An interesting observation of this study

382

was the increased level of Fis1, which is known to enable mitochondrial fission, as well as the

383

decreased levels of Mfn1/2 and Opa1, which facilitate mitochondrial fusion. Although the role of

384

mitochondrial fission and fusion during the period up to the implantation stage has not been

385

well-defined, mouse embryonic fibroblasts (MEFs), as well as ES and TS cells from mouse knockout

386

models displaying deficient fission and fusion, showed evident aberrations in mitochondrial mass and

387

morphology (elongated or fragmented mitochondria) [37-39]. Thus, we predict that the resumption of

388

mitochondrial growth by the implantation stage is not only facilitated by increased mitochondrial

389

biogenesis and mtDNA replication, but also by the appropriate balance between fission and fusion (Fig.

390

7A). In addition, deficient fission and fusion result in disrupted mitochondrial functions, distribution,

391

and apoptosis, with deficient embryos dying in midgestation [37-39]. These facts imply that fission

392

and fusion dynamics may also be involved in many aspects of mitochondrial modifications during this

393

critical period.

394

Enhanced mitochondrial functions not only provide more ATP but also produce more ROS.

395

Energy metabolism in mitochondria, including TCA, β-oxidation, and (especially) OXPHOS, is the

396

major source of ROS. Our results regarding the GSH/GSSG ratio suggested that embryos face

397

increased oxidative stress soon after implantation. Mitochondria are both a major source and a target

398

of ROS [40]. Mitochondrial glutathione (mGSH), which is synthesized in the cytosol and then

399

imported into the mitochondrial matrix, is thought to be the key survival antioxidant for maintaining

400

the appropriate mitochondrial redox environment [41]. Interestingly, the temporal pattern of GSH

401

biosynthesis is similar to the dynamics of the resumption of mitochondrial biogenesis by the

402

implantation stage [42]. Thus, increased GSH biosynthesis may be a physiological response to

403

enhanced mitochondrial functions and the consequent increase in ROS levels. Indeed, we observed an

404

evident activation of many key MtGs involved in GSH biosynthesis and metabolism. The expression

405

of Gclc and Gs (Glul), which code for the enzymes that catalyze the first (rate-limiting) and second

406

steps of glutathione biosynthesis, increased during implantation. We also found that Nrf1 and Nrf2

407

were upregulated; these well-known transcription factors can activate GSH-based antioxidative genes,

408

including Gclc, Gss, Glrxs, and Gpxs [43, 44]. Consistently, Keap1 and Cul3, which code for proteins

409

that inactivate Nrf2 by forming a KAP1/CUL3/NRF2 complex in the cytosol, were downregulated

410

during the transition from E3.5 blastocysts to E7.5 epiblasts. Previous studies in normal somatic cells

411

and cancer cells indicated that the inhibition of either Keap1 or Cul3 increases the nuclear

412

accumulation of Nrf2, leading to the activation of Nrf2-dependent antioxidative gene expression [45,

413

46]. Although the role of ROS in uncoupling NRF2 from the KAP1/CUL3 complex and thereby

414

upregulating the transcription of GSH-based antioxidative genes has been well documented in many

415

somatic cell types [47], the involvement of ROS in the Nrf2-induced activation of the GSH system has

416

never been demonstrated in early embryogenesis. Our data suggest that the enhancement of

417

mitochondrial functions in embryos by the implantation stage may increase the production of cellular

418

ROS, activating the expression of the GSH-based antioxidative system via the Nrf2-Keap1 signaling

419

pathway to protect mitochondria against oxidative stress (Fig. 7B). The transition from E3.5

420

blastocysts to E7.5 epiblasts was also characterized by a marked increase in the expression of many

421

Gpxs, which are the genes coding for proteins that catalyze the reduction of hydrogen peroxide and

422

lipid peroxides by GSH. Among these, Gpx3 and Gpx4 are well-known mitochondrial Gpx genes that

423

are essential for oxidative homeostasis in mitochondria [48, 49].

424

We also found many apoptotic genes undergoing significant changes during early development.

425

TUNEL analysis showed that E7.5 epiblasts displayed a high rate of apoptosis compared to E3.5

426

blastocysts, consistently to previously reported observations in hESCs undergoing differentiation [50].

427

Compared with MtGs associated with ‘Mitochondrial biogenesis and mtDNA replication’ and

428

‘GSH-based antioxidative response’, the gene expression dynamics of ‘Apoptotic mitochondrial

429

change’ showed a complicated pattern. Indeed, properly controlled mitochondrial apoptotic signaling

430

is not harmful to early embryos; in contrast, it is considered significant for the survival and

431

differentiation of early embryos [50, 51]. A PPI analysis showed that the term ‘Apoptotic

432

mitochondrial change’ was highly associated with ‘GSH-based antioxidative response’ and

433

‘Mitochondrial biogenesis’. Similarly, REVIGO analysis showed that ‘Apoptotic mitochondrial

434

change’ interacted significantly with ‘Oxoacid metabolism’. These observations imply that the

435

apoptotic and survival signals may be balanced by the GSH-based stress response (Fig. 7B).

436

Our data suggest that, in addition to successful implantation, mitochondrial functions or dynamics

437

may be also involved in early organogenesis. Using the MGI database, we found that many of the

438

MtGs whose expression was significantly altered (P < 0.05, FC >2) were associated with embryonic

439

death or aberrant development. At E7.5, when embryos are in an advanced stage of the gastrulation

440

process and have a formed neural plate, many of these MtGs were annotated with phenotypes such as

441

‘Decreased neuron number’, ‘Abnormal neural crest cell apoptosis’, and ‘Abnormal neural tube

442

morphology/development’. From E7.5 to E10.5, the neural plate folds to form the neural tube and the

443

brain, while structures and organs such as somites, the heart, and the limb buds start developing. MtGs

444

at this stage were annotated with phenotypes that included ‘Abnormal nervous system morphology’,

445

‘Abnormal eye development’, ‘Abnormal heart morphology’, and ‘Abnormal limb morphology’. It is

446

noteworthy that the MAPK and Wnt signaling pathways were tightly clustered with MtGs that are

447

responsible for mitochondrial functions. Previous studies in mouse and human ES cells indicated that

448

mitochondrial apoptotic signals might contribute to differentiation programs, probably via the MAPK

449

and Wnt signaling pathways [50, 51].

450

Our next step was to investigate whether the well-orchestrated mitochondrial functions and

451

dynamics taking place by the implantation stage are regulated by epigenetic modifications. Coinciding

452

with the dramatic mitochondrial changes during early development, embryos undergo de novo DNA

453

methylation. Aberrant de novo DNA methylation caused by deleting or mutating Dnmts in mice or ES

454

cells results in embryonic lethality and failed differentiation [52, 53]. Thus, we decided to investigate

455

the possible functional association of these two early developmental events. Our data showed that

456

most MtGs (clusters 1 and 2) underwent de novo methylation in their promoter regions during early

457

development. These observations were very similar to those reported by Borgel et al. [19] and Smith et

458

al. [54], who demonstrated that embryos undergo global promoter de novo methylation in the period

459

between the preimplantation and the postimplantation stages. In our study, the integrated analysis of

460

MtGs that underwent promoter de novo methylation revealed that approximately half of them

461

displayed reduced expression levels. Many of these MtGs are significantly associated with

462

mitochondrial functions and diseases. These results demonstrate that de novo DNA methylation, one of

463

the most important epigenetic events during mammalian embryogenesis, may have considerable

464

functional consequences for the mitochondrial modifications that are essential for embryonic survival

465

and growth.

466

Another finding was that the mtDNA of both E7.5 and E10.5 embryos showed significant

467

enrichment of MeDIP reads. This novel finding suggests that mtDNA methylation, which has been

468

reported in many somatic cell types [20, 55, 56], can be observed as early as immediately after

469

implantation. We also showed that this enrichment in early embryos was not limited to regions that

470

coded for enzymes essential for OXPHOS (ND5, ND6), but also in noncoding regions (12S rRNA) and

471

a control region (D-loop). In addition, in our study mtDNA methylation appeared to be independent of

472

CpG density, as it was detected in both high- and low-CpG regions. This result is consistent with the

473

findings of Bellizzi et al., who reported that both CpG and non-CpG regions of mtDNA can be

474

methylated in human blood and cultured cells [56]. To the best of our knowledge, this is the first study

475

to show that mtDNA becomes methylated in early embryos and that mt-DNMT1 can translocate into

476

the mitochondria matrix and bind to mtDNA for maintaining mtDNA methylation by the implantation

477

stage. Indeed, our unpublished data indicate that mtDNA undergoes a process similar to the de novo

478

DNA methylation of the nuclear genome, thereby demonstrating that highly dynamic DNA

479

methylation modifications are critical not only for the nuclear genome but also for the mitochondrial

480

genome.

481

Overall, our study provided high-resolution dynamic expression profiles of MtGs during early

482

embryogenesis. Our data indicate that orchestrated MtG expression is critical for the ‘embryonic shift’

483

in energy metabolism and the resumption of mitochondrial biogenesis, which is of prime importance

484

for successful embryonic implantation and survival. In addition, the mitochondria-mediated balance

485

between apoptosis and the stress response plays a significant role in the differentiation and

486

organogenesis of postimplantation embryos. Our data also suggest that promoter de novo methylation

487

of nDNA-encoded MtGs may have considerable functional consequences for mitochondrial

488

modifications. Finally, mtDNA methylation in somatic cells can be observed as early as immediately

489

after implantation. Our study can serve as a comprehensive reference for the further study of

490

mitochondrial modulations during early mammalian embryogenesis and provides a new insight into

491

the epigenetic mechanisms modulating mitochondrial functions and dynamics at this stage.

492

Acknowledgments

493

This work was supported by grants from the National Natural Science Foundation of China (No.

494

31672426 and 31472092) and the National High-Tech R&D Program (2013AA102506).

495 496

References:

497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540

1. Harvey A, Gibson T, Lonergan T, Brenner C. Dynamic regulation of mitochondrial function in preimplantation embryos and embryonic stem cells. Mitochondrion 2011; 11: 829-838. 2. Ramalho-Santos J, Varum S, Amaral S, Mota PC, Sousa AP, Amaral A. Mitochondrial functionality in reproduction: from gonads and gametes to embryos and embryonic stem cells. Hum Reprod Update 2009; 15: 553-572. 3. Mootha VK, Bunkenborg J, Olsen JV, Hjerrild M, Wisniewski JR, Stahl E, Bolouri MS, Ray HN, Sihag S, Kamal M, Patterson N, Lander ES, Mann M. Integrated analysis of protein composition, tissue diversity, and gene regulation in mouse mitochondria. Cell 2003; 115: 629-640. 4. Jian B, Yang S, Chen D, Zou L, Chatham JC, Chaudry I, Raju R. Aging influences cardiac mitochondrial gene expression and cardiovascular function following hemorrhage injury. Mol Med 2011; 17: 542-549. 5. Hsiao CP, Wang D, Kaushal A, Saligan L. Mitochondria-related gene expression changes are associated with fatigue in patients with nonmetastatic prostate cancer receiving external beam radiation therapy. Cancer Nurs 2013; 36: 189-197. 6. Ren L, Wang Z, An L, Zhang Z, Tan K, Miao K, Tao L, Cheng L, Zhang Z, Yang M, Wu Z, Tian J. Dynamic comparisons of high-resolution expression profiles highlighting mitochondria-related genes between in vivo and in vitro fertilized early mouse embryos. Hum Reprod 2015; 30: 2892-2911. 7. Takasugi M, Yagi S, Hirabayashi K, Shiota K. DNA methylation status of nuclear-encoded mitochondrial genes underlies the tissue-dependent mitochondrial functions. BMC Genomics 2010; 11: 481. 8. Kelly RD, Mahmud A, McKenzie M, Trounce IA, St John JC. Mitochondrial DNA copy number is regulated in a tissue specific manner by DNA methylation of the nuclear-encoded DNA polymerase gamma A. Nucleic Acids Res 2012; 40: 10124-10138. 9. Hashemi M, Rezaei H, Kaykhaei MA, Taheri M. Association of the Promoter Methylation of Mitochondrial Transcription Factor A With Susceptibility to Metabolic Syndrome. Gene Cell Tissue 2014; 1: e18308. 10. Gemma C, Sookoian S, Dieuzeide G, Garcia SI, Gianotti TF, Gonzalez CD, Pirola CJ. Methylation of TFAM gene promoter in peripheral white blood cells is associated with insulin resistance in adolescents. Mol Genet Metab 2010; 100: 83-87. 11. Barres R, Osler ME, Yan J, Rune A, Fritz T, Caidahl K, Krook A, Zierath JR. Non-CpG methylation of the PGC-1alpha promoter through DNMT3B controls mitochondrial density. Cell Metab 2009; 10: 189-198. 12. Downs KM, Davies T. Staging of gastrulating mouse embryos by morphological landmarks in the dissecting microscope. Development 1993; 118: 1255-1266. 13. Sui L, An L, Tan K, Wang Z, Wang S, Miao K, Ren L, Tao L, He S, Yu Y, Nie J, Liu Q, Xing L, Wu Z, Hou Z, Tian J. Dynamic proteomic profiles of in vivo- and in vitro-produced mouse postimplantation extraembryonic tissues and placentas. Biol Reprod 2014; 91: 155. 14. Giritharan G, Talbi S, Donjacour A, Di Sebastiano F, Dobson AT, Rinaudo PF. Effect of in vitro fertilization on gene expression and development of mouse preimplantation embryos. REPRODUCTION 2007; 134: 63-72. 15. Huang W, Khatib H. Comparison of transcriptomic landscapes of bovine embryos using RNA-Seq. BMC Genomics 2010; 11: 711. 16. Kepkova KV, Vodicka P, Toralova T, Lopatarova M, Cech S, Dolezel R, Havlicek V, Besenfelder U, Kuzmany A, Sirard MA, Laurincik J, Kanka J. Transcriptomic analysis of in vivo and in vitro produced bovine embryos revealed

a developmental change in cullin 1 expression during maternal-to-embryonic transition.

Theriogenology 2011; 75: 1582-1595.

541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585

17. Nie J, An L, Miao K, Hou Z, Yu Y, Tan K, Sui L, He S, Liu Q, Lei X, Wu Z, Tian J. Comparative analysis of dynamic proteomic profiles between in vivo and in vitro produced mouse embryos during postimplantation period. J Proteome Res 2013; 12: 3843-3856. 18. Zitka O, Skalickova S, Gumulec J, Masarik M, Adam V, Hubalek J, Trnkova L, Kruseova J, Eckschlager T, Kizek R. Redox status expressed as GSH:GSSG ratio as a marker for oxidative stress in paediatric tumour patients. Oncol Lett 2012; 4: 1247-1253. 19. Borgel J, Guibert S, Li Y, Chiba H, Schubeler D, Sasaki H, Forne T, Weber M. Targets and dynamics of promoter DNA methylation during early mouse development. Nat Genet 2010; 42: 1093-1100. 20. Shock LS, Thakkar PV, Peterson EJ, Moran RG, Taylor SM. DNA methyltransferase 1, cytosine methylation, and cytosine hydroxymethylation in mammalian mitochondria. Proc Natl Acad Sci U S A 2011; 108: 3630-3635. 21. Geng X, Oliver G. Elucidating the molecular characteristics of organogenesis in human embryos. Genome Biol 2010; 11: 130. 22. Jornayvaz FR, Shulman GI. Regulation of mitochondrial biogenesis. Essays Biochem 2010; 47: 69-84. 23. Lee HC, Wei YH. Mitochondrial biogenesis and mitochondrial DNA maintenance of mammalian cells under oxidative stress. Int J Biochem Cell Biol 2005; 37: 822-834. 24. Hance N, Ekstrand MI, Trifunovic A. Mitochondrial DNA polymerase gamma is essential for mammalian embryogenesis. Hum Mol Genet 2005; 14: 1775-1783. 25. Larsson NG, Wang J, Wilhelmsson H, Oldfors A, Rustin P, Lewandoski M, Barsh GS, Clayton DA. Mitochondrial transcription factor A is necessary for mtDNA maintenance and embryogenesis in mice. Nat Genet 1998; 18: 231-236. 26. Humble MM, Young MJ, Foley JF, Pandiri AR, Travlos GS, Copeland WC. Polg2 is essential for mammalian embryogenesis and is required for mtDNA maintenance. Hum Mol Genet 2013; 22: 1017-1025. 27. Kazak L, Reyes A, Holt IJ. Minimizing the damage: repair pathways keep mitochondrial DNA intact. Nat Rev Mol Cell Biol 2012; 13: 659-671. 28. He J, Cooper HM, Reyes A, Di Re M, Sembongi H, Litwin TR, Gao J, Neuman KC, Fearnley IM, Spinazzola A, Walker JE, Holt IJ. Mitochondrial nucleoid interacting proteins support mitochondrial protein synthesis. Nucleic Acids Res 2012; 40: 6109-6121. 29. Arakaki N, Nishihama T, Kohda A, Owaki H, Kuramoto Y, Abe R, Kita T, Suenaga M, Himeda T, Kuwajima M, Shibata H, Higuti T. Regulation of mitochondrial morphology and cell survival by Mitogenin I and mitochondrial single-stranded DNA binding protein. Biochim Biophys Acta 2006; 1760: 1364-1372. 30. Christie DA, Lemke CD, Elias IM, Chau LA, Kirchhof MG, Li B, Ball EH, Dunn SD, Hatch GM, Madrenas J. Stomatin-like protein 2 binds cardiolipin and regulates mitochondrial biogenesis

and function.

Mol Cell Biol 2011; 31: 3845-3856. 31. Rampazzo C, Fabris S, Franzolin E, Crovatto K, Frangini M, Bianchi V. Mitochondrial thymidine kinase and the enzymatic network regulating thymidine triphosphate pools in cultured human cells. J Biol Chem 2007; 282: 34758-34769. 32. Artuso L, Romano A, Verri T, Domenichini A, Argenton F, Santorelli FM, Petruzzella V. Mitochondrial DNA metabolism in early development of zebrafish (Danio rerio). Biochim Biophys Acta 2012; 1817: 1002-1011. 33. Bossy-Wetzel E, Barsoum MJ, Godzik A, Schwarzenbacher R, Lipton SA. Mitochondrial fission in apoptosis, neurodegeneration and aging. Curr Opin Cell Biol 2003; 15: 706-716. 34. Santel A, Fuller MT. Control of mitochondrial morphology by a human mitofusin. J Cell Sci 2001; 114: 867-874. 35. Chen Y, Liu Y, Dorn GN. Mitochondrial fusion is essential for organelle function and cardiac homeostasis.

586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630

Circ Res 2011; 109: 1327-1331. 36. Ishihara N, Fujita Y, Oka T, Mihara K. Regulation of mitochondrial morphology through proteolytic cleavage of OPA1. EMBO J 2006; 25: 2966-2977. 37. Chen H, Detmer SA, Ewald AJ, Griffin EE, Fraser SE, Chan DC. Mitofusins Mfn1 and Mfn2 coordinately regulate mitochondrial fusion and are essential for embryonic development. J Cell Biol 2003; 160: 189-200. 38. Chen H, Chomyn A, Chan DC. Disruption of fusion results in mitochondrial heterogeneity and dysfunction. J Biol Chem 2005; 280: 26185-26192. 39. Ishihara N, Nomura M, Jofuku A, Kato H, Suzuki SO, Masuda K, Otera H, Nakanishi Y, Nonaka I, Goto Y, Taguchi N, Morinaga H, Maeda M, Takayanagi R, Yokota S, Mihara K. Mitochondrial fission factor Drp1 is essential for embryonic development and synapse formation in mice. Nat Cell Biol 2009; 11: 958-966. 40. Reinecke F, Smeitink JA, van der Westhuizen FH. OXPHOS gene expression and control in mitochondrial disorders. Biochim Biophys Acta 2009; 1792: 1113-1121. 41. Mari M, Morales A, Colell A, Garcia-Ruiz C, Fernandez-Checa JC. Mitochondrial glutathione, a key survival antioxidant. Antioxid Redox Signal 2009; 11: 2685-2700. 42. Hansen JM, Harris C. Glutathione during embryonic development. Biochim Biophys Acta 2015; 1850: 1527-1542. 43. Yang H, Magilnick N, Lee C, Kalmaz D, Ou X, Chan JY, Lu SC. Nrf1 and Nrf2 regulate rat glutamate-cysteine ligase catalytic subunit transcription indirectly via NF-kappaB and AP-1. Mol Cell Biol 2005; 25: 5933-5946. 44. Chepelev NL, Zhang H, Liu H, McBride S, Seal AJ, Morgan TE, Finch CE, Willmore WG, Davies KJ, Forman HJ. Competition of nuclear factor-erythroid 2 factors related transcription factor isoforms, Nrf1 and Nrf2, in antioxidant enzyme induction. Redox Biol 2013; 1: 183-189. 45. Cullinan SB, Gordan JD, Jin J, Harper JW, Diehl JA. The Keap1-BTB protein is an adaptor that bridges Nrf2 to a Cul3-based E3 ligase: oxidative stress sensing by a Cul3-Keap1 ligase. Mol Cell Biol 2004; 24: 8477-8486. 46. Wang R, An J, Ji F, Jiao H, Sun H, Zhou D. Hypermethylation of the Keap1 gene in human lung cancer cell lines and lung cancer tissues. Biochem Biophys Res Commun 2008; 373: 151-154. 47. Hayes JD, Dinkova-Kostova AT. The Nrf2 regulatory network provides an interface between redox and intermediary metabolism. Trends Biochem Sci 2014; 39: 199-218. 48. Passaia G, Fonini LS, Caverzan A, Jardim-Messeder D, Christoff AP, Gaeta ML, de Araujo MJ, Margis R, Margis-Pinheiro M. The mitochondrial glutathione peroxidase GPX3 is essential for H2O2 homeostasis and root and shoot development in rice. Plant Sci 2013; 208: 93-101. 49. Savaskan NE, Ufer C, Kuhn H, Borchert A. Molecular biology of glutathione peroxidase 4: from genomic structure to developmental expression and neural function. Biol Chem 2007; 388: 1007-1017. 50. Qin H, Yu T, Qing T, Liu Y, Zhao Y, Cai J, Li J, Song Z, Qu X, Zhou P, Wu J, Ding M, Deng H. Regulation of apoptosis and differentiation by p53 in human embryonic stem cells. J Biol Chem 2007; 282: 5842-5852. 51. Duval D, Trouillas M, Thibault C, Dembele D, Diemunsch F, Reinhardt B, Mertz AL, Dierich A, Boeuf H. Apoptosis and differentiation commitment: novel insights revealed by gene profiling studies in mouse embryonic stem cells. Cell Death Differ 2006; 13: 564-575. 52. Okano M, Bell DW, Haber DA, Li E. DNA methyltransferases Dnmt3a and Dnmt3b are essential for de novo methylation and mammalian development. Cell 1999; 99: 247-257. 53. Liao J, Karnik R, Gu H, Ziller MJ, Clement K, Tsankov AM, Akopian V, Gifford CA, Donaghey J, Galonska C, Pop R, Reyon D, Tsai SQ, Mallard W, Joung JK, Rinn JL, Gnirke A, Meissner A. Targeted disruption of DNMT1, DNMT3A and DNMT3B in human embryonic stem cells. Nat Genet 2015; 47: 469-478.

631 632 633 634 635 636

54. Smith ZD, Chan MM, Mikkelsen TS, Gu H, Gnirke A, Regev A, Meissner A. A unique regulatory phase of DNA methylation in the early mammalian embryo. Nature 2012; 484: 339-344. 55. Dzitoyeva S, Chen H, Manev H. Effect of aging on 5-hydroxymethylcytosine in brain mitochondria. Neurobiol Aging 2012; 33: 2881-2891. 56. Bellizzi D, D'Aquila P, Scafone T, Giordano M, Riso V, Riccio A, Passarino G. The control region of mitochondrial DNA shows an unusual CpG and non-CpG methylation pattern. DNA Res 2013; 20: 537-547.

637 638

639

Figure legends

640

Figure 1. General dynamic profiling of mitochondria-related genes (MtGs) during early development.

641

(A) Hierarchical clustering analysis based on the dynamic expression levels of MtGs of embryos from

642

pre- to postimplantation. Normalized expression levels (RPKM) are represented by different colors:

643

red indicates high abundance and green indicates low abundance. The number of MtGs in each cluster

644

is presented on the right side. General trends in expression changes are indicated by the average

645

RPKM value of MtGs in each cluster. (B) Fold change (FC) distribution analysis of MtGs between

646

E3.5 blastocysts to E7.5 epiblasts as well as between E7.5 epiblasts and E10.5 embryos. Only MtGs

647

with statistically significant changes in expression were included (P < 0.05). (C) Venn diagrams of

648

upregulated (left panel, red) and downregulated (right panel, green) MtGs for two developmental

649

transitions. Tables show the major functions of MtGs that are consistently upregulated (left) or

650

downregulated (right) during early development.

651

Figure 2. General functional profiling of MtGs whose expression significantly changed (P < 0.05,

652

FC >2) during early development. (A-B) Classification of GO terms based on the functional

653

annotation of biological processes (BPs) enriched in the transition from E3.5 blastocysts to E7.5

654

epiblasts, as well as from E7.5 epiblasts to E10.5 embryos. The left ordinate represents the number of

655

enriched MtGs corresponding to each term and the right ordinate represents the enrichment score

656

(defined as −Log10 P-value). (C-D) Graph visualization of enriched BPs in A and B based on REVIGO

657

analysis. Bubble color indicates P-value. Functionally associated BPs are linked. (E-F) Interaction

658

networks of MtGs during early development created by a web-based search of the STRING database.

659

Boxed regions represent tightly interconnected functional clusters.

660

Figure 3. Mitochondrial modifications during implantation stage. (A) Model illustrating the

661

‘embryonic shift’ from anaerobic to aerobic metabolism, as well as the resumption of mitochondrial

662

biogenesis and the enhanced response to oxidative stress. Red text indicates upregulated MtGs, while

663

green text indicates downregulated MtGs. (B) Fold changes of representative MtGs involved in

664

essential mitochondrial modifications by the implantation stage. (C) Relative expression levels of

665

some representative MtGs in E3.5 blastocysts and E7.5 epiblasts determined by qRT-PCR. (D)

666

Interaction networks of MtGs responsible for essential mitochondrial modifications taking place

667

during the implantation stage. (E) Boxplot quantitative comparison of mtDNA copy numbers between

668

E10.5 embryos with higher and lower survival rates. (F) Ratios of reduced glutathione (GSH) to

669

oxidized glutathione (GSSG) in E3.5 blastocysts and E7.5 epiblasts; **P < 0.01. (G) Representative

670

images of apoptosis detected by TUNEL in E3.5 blastocysts and E7.5 epiblasts. Rightmost panels:

671

higher magnification of boxed regions in E7.5 epiblasts. The nuclei of blastomeres were labeled with

672

DAPI. Representative TUNEL-positive apoptotic nuclei are indicated with arrows. Scale bars, 100

673

μm.

674

Figure 4. Schematic diagram of the functional associations among MtGs that are annotated with

675

aberrant organogenesis during the postimplantation stage. Colored boxes indicate mitochondrial

676

modifications that may be essential for embryonic survival and normal organogenesis during the

677

implantation and postimplantation periods.

678

Figure 5. Dynamics of promoter DNA methylation of nDNA-encoded MtGs and functional

679

consequences during early development (A) Hierarchical clustering analysis based on dynamic

680

methylation levels of MtGs during the transition from pre- to postimplantation embryos. Normalized

681

methylation levels are represented by different colors: red indicates relatively hypermethylated

682

promoters and black indicates relatively hypomethylated promoters. The numbers of MtGs in each

683

cluster are displayed on the right side of the plot. Graphs next to the plot display promoter methylation

684

levels during early development. The red line (upper graph) indicates total normalized methylation

685

values while blue lines (three lower graphs) correspond to average normalized methylation levels in

686

each cluster. (B) Scatter plots of expression changes versus mean differences of methylation for the

687

MtGs of cluster 1 (upper plot, boxed by green dotted line) and cluster 2 (lower plot, boxed by blue

688

dotted line) whose promoters underwent de novo methylation before E7.5 and E10.5, respectively.

689

Green dots indicate MtGs that are downregulated during the process of promoter de novo methylation,

690

while red dots indicate upregulated MtGs (C-D) Functional profiling of MtGs corresponding to the

691

green dots of the shaded regions of the scatter plots in (B) and may be negatively regulated by

692

promoter de novo methylation taking place between E3.5 and E7.5 or between E7.5 and E10.5.

693

Figure 6. MtDNA methylation in postimplantation embryos. (A) Expression dynamics of Dnmt1

694

during early development as assessed by RNA-seq. (B) Relative expression levels of total Dnmt1 and

695

mt-Dnmt1 during early development as determined by qRT-PCR. (C) MtDNA methylation levels are

696

presented as enriched MeDIP-seq reads (y-axis) on mtDNA positions (x-axis) in E7.5 epiblasts (upper

697

plot) and E10.5 embryos (lower plot). (D) Relative methylation levels at selected regions of mtDNA as

698

detected by MeDIP-qPCR in E7.5 epiblasts (upper plot) and E10.5 embryos (lower plot). (E) Relative

699

enrichment of DNMT1 at selected regions (Nd5 and Nd6) of mtDNA as determined by CHIP-qPCR in

700

E6.5 embryos. Values indicated by different letters are significantly different (P < 0.05). *P < 0.05;

701

**P < 0.01.

702

Figure 7. Schematic diagrams illustrating essential mitochondrial modifications during early

703

development. (A) The resumption of mitochondrial biogenesis by the implantation stage appears to be

704

facilitated not only by increased mitochondrial biogenesis and mtDNA replication but also by the

705

skewed balance between fission and fusion. (B) Enhanced aerobic metabolism by the implantation

706

stage, especially OXPHOS, leads to increased ROS production, which appears be the inducer of GSH

707

biosynthesis. The apoptotic and survival signals in early embryos may be balanced by the GSH-based

708

stress response.

709

710

Table 1. Enriched KEGG pathways (Top 15) for MtGs significantly changed during early

711

development E7.5 epiblasts vs. E3.5 blastocysts

KEGG terms

E10.5 embryos vs. E7.5 epiblasts

Count

P-value

KEGG terms

Count

P-value

22

6.56E-16

Alzheimer's disease

24

2.45E-09

Fatty acid metabolism

17

1.60E-11

14

2.57E-09

Alzheimer's disease

31

8.82E-11

Huntington's disease

24

2.74E-09

Butanoate metabolism

14

1.63E-09

Parkinson's disease

18

2.92E-07

15

3.27E-09

Oxidative phosphorylation

17

1.12E-06


28

1.03E-08

Propanoate metabolism

8

2.03E-05

Propanoate metabolism

12

1.83E-08

Butanoate metabolism

8

8.60E-05

Pyruvate metabolism

13

7.36E-08

9

2.41E-04

PPAR signaling pathway

17

1.45E-07

ABC transporters

8

3.11E-04

Parkinson's disease

22

1.51E-07

Fatty acid metabolism

8

3.11E-04

Citrate cycle (TCA cycle)

11

3.37E-07

8

3.57E-04

14

5.22E-07

10

0.0010

15

7.23E-07

7

0.0019

Tryptophan metabolism

10

3.48E-05

9

0.0022

beta-Alanine metabolism

7

2.39E-04

6

0.0023

Arginine and proline metabolism

Valine, leucine and isoleucine degradation

Amyotrophic lateral sclerosis (ALS) Glycolysis / Gluconeogenesis

712

Arginine and proline metabolism

Amyotrophic lateral sclerosis (ALS)

Valine, leucine and isoleucine degradation Apoptosis Steroid hormone biosynthesis PPAR signaling pathway Glycine, serine and threonine metabolism

92

Cluster 5

203

Cluster 6

50 40 30 20 10

Cluster 5

60 45 30 15 0

Log 2 Fold change

Group 1 Cluster 6

Group 3

144

Cluster 4

Group 3

Cluster 4

10 38

6

107 116 223 248 90 74

2 -2

261 920 188

24

-6

-10 10

Expression levels (RPKM)

Cluster 3

Group 2

100 80 60 40 20

121

172

Cluster 2

Cluster 3

120 100 80 60 40

B

Cluster 2 100 80 60 40 20

Log 2 Fold change

Cluster 1

Group 1

Cluster 1 200 170 140 110 80 50

Group 2

E10.5

306

E7.5


E3.5


A

21

6

175 54 100 260 281 74 42

2 -2 -6

843 127

11

-10 -10 -5 0 5 10 15 Log2 average RPKM

C

Functions

Consistently upregulated MtGs (FC>2)

Apoptosis/ proliferation/ differentiation

Nfatc4, Prkca, Smad3, Pak7, Ndn, Ak1, Casp8, Cyp1b1, Stmn1, Dpysl2, Bcl2, Wnt1, Shh, Cav1, Nos3, Mtch1, Stard13

Lipid metabolism

Pla2g4a, Abcd2, Mlycd, Hmgcs2

E7.5 vs. E3.5 193 58

151 31

E10.5 vs. E7.5

94 Mpo, Glrx2, Idh1, Mgst1

transmembra ne transport

Scn1b, Slc8a2, Slc25a25, Cacna1b

Consistently down-regulated mito-genes (FC>2)

Apoptosis/ proliferation/ differentiation

Cyct, Dmc1, Lgals12, Mtus1, Spns1, Ndufs1, Spast

Genetic information

Dmc1, Ppargc1b

Mito-membrane organization

Mtx1, Dsp, Immt, Slc25a41, Efhd1, Timm9, Cacna1d, Myh9

Glucose/ lipid/ amino acis metabolism

Echdc2, Atp5g1, Cox7a2l, Dmgdh, Atg7, Ass1

E7.5 vs. E3.5

111 Anti-oxidation

Functions

E10.5 vs. E7.5

Figure 1

No. of enriched genes

35 30 25 20 15 10 5 0

60 40 20 0

C

D

E

F

Mitochondrial translation

β-oxidation & TCA

Enrichment score


80

Enrichment score


B 50

25

40

20

30

15

20

10

10

5

0

0

Enrichment score

A

OXPHOS & mitochondria translation & transmembrane transport

Electron transfer chain

Amino acid metabolism Transmembrane transport

OXPHOS

Response to ROS Apoptosis/ proliferation

TCA & β-oxidation & Amino acid metabolism

Mitochondrial biogenesis

Figure 2

Response to ROS & apoptosis & proliferation & differentiation Prostaglandin metabolism

E3.5

B

E7.5

Hk2, Pdha1, Pdha2

Glycolysis Increased ROS OXPHOS

Atp5a1, Atp5b, Atp5d, Atp5e, Atp5f1, Atp5g1, Atp5g2, Atp5h, Atp5j2, Atp5s

Log2 Fold change form E3.5 to E7.5

A

Response to oxidative stress

Gclc, Gs, Gpx3, Gpx7, Glrx5, Gstm1, Slc25a11, Cul3, Keap1



Polg, Slc25a33, Polg2, Tfam, Ssbp1, Dhodh, Fis1, Mfn1, Mfn2

5

MtGs repressing mitochondrial biogenesis

8 7 6 5 4 3 2 1 0 -1 -2 -3

GSH-based antioxidative response

MtGs repressing GSH synthesis

E3.5

4

E7.5

3


Relative expression levels to E3.5

Mitochondrial biogenesis & mtDNA replication

Cpt1b, Cpt1c, Cpt2

β-oxidation

C

6 5 4 3 2 1 0 -1 -2 -3

2 1 0

8 7 6 5 4 3 2 1 0 -1 -2

Apoptotic mitochondrial change

GSH-based antioxidative response

D

E Apoptotic mitochondrial change


mtDNA copy number

1500 1200 900 600 300 0

80

G

E7.5 epiblast

Boxed region a

DAPI

100 60 40 20 0 E3.5 (n = 30)

E7.5 (n = 20)

DAPI/ TUNEL

F

Ratio of GSH/ GSSG

E3.5 blastocysts

**

Figure 3

a

b

b

Hk2: abnormal

Lrp5: abnormal eye

neuronal precursor proliferation

development; abnormal limb morphology

Cybb:

Mthfd1: abnormal

cardiac hypertrophy

Mthfd1l:

abnormal neural tube morphology

heart morphology; abnormal neural tube closure

Glucose/ lipid/ amino acid metabolism

vascular endothelial cell morphology

Dnaja3:

enlarged heart; cardiac interstitial fibrosis

Nos3: abnormal heart

morphology; abnormal limb morphology; abnormal neuron morphology; abnormal vascular branching morphogenesis

Cav1: abnormal astrocyte morphology; cardiomyopathy

Drd2: abnormal neuronal

Mitochondrial biogenesis & mtDNA replication

morphology; abnormal sensory capabilities

Shh : abnormal eye Polg: dilated

Abcd2 : abnormal nervous

heart left ventricle

system physiology; abnormal nervous system physiology

development; abnormal heart morphology; abnormal limb morphology; abnormal nervous system morphology

Hspa5 : abnormal

Slc8a2 :

abnormal neuron physiology

bone morphology

Apoptotic changes

abnormal sensory capabilities

abnormal neuron physiology

neuron differentiation/ proliferation

Atp7a: abnormal limb

Transmembrane transport

Cacna1b :

Scn2a1:

migration; abnormal neuron physiology

Nos2: abnormal heart

Abcd1: abnormal lens morphology; abnormal nervous system physiology

Ptrf: abnormal

Pink1: abnormal

neuron physiology; abnormal sensory capabilities

Ucp2: decreased

Vdr : abnormal heart ventricle morphology; abnormal limb bone morphology

Rnf7 : abnormal neuronal

neuron apoptosis

Casp2: abnormal

migration

neuron apoptosis

Response to stress

Bbc3 : abnormal heart left Glrx2 : abnormal

ventricle pressure; abnormal neuron physiology

eye physiology

Sod1 : neurodegeneration

Fxn : neurodegeneration Gstm1 : increased neuron apoptosis

Apaf1 : abnormal neural tube

morphology; abnormal neuronal migration/differentation; increased neuron apoptosis

Figure 4

Log2 Fold change (RNA-seq)

500

Cluster 1

0

Cluster 2 18 15 12 9 6 3 0

10 8 6 4 2

D

0

4

4.5

5

-3

48

Cluster 2

6

102

3 0

4

4.5

5

5.5

6

-3

85

-6

Log2 Fold change (MeDIP-seq)

Cluster 1 10

8

8

6

6

4

4 2

2

0

0

-Log 10 P value

4.0 3.0 2.0 1.0 0.0

0

Log2 Fold change (MeDIP-seq)

Log2 Fold change (RNA-seq)

3.5 2.8 2.1 1.4 0.7 0.0

40

3

-6

Cluster 2

3.0 2.4 1.8 1.2 0.6 0.0

-Log 10 P value

C


Cluster 3

Cluster 1

6

1000

Cluster 3

69

212

Cluster 2

B

Total 1500


Cluster 1

Ave. normalized Ave. normalized Ave. normalized methylation levels methylation levels methylation levels

103

Total normalized methylation levels

A

KEGG pathway (Cluster 2)

Count

P-value

KEGG pathway (Cluster 1)

Count

Pvalue

Aminoacyl-tRNA biosynthesis

5

0.0002

Glyoxylate and dicarboxylate metabolism

3

0.001

Alzheimer's disease

7

0.002

Alzheimer's disease

5

0.002

Arginine & proline metabolism

4

0.008


5

0.002


6

0.012

Parkinson's disease

4

0.008

Alanine, aspartate ,glutamate metabolism

3

0.022

Amyotrophic lateral sclerosis

3

0.014

Phenylalanine, tyrosine and tryptophan biosynthesis

2

0.038

Apoptosis

3

0.032

Figure 5

0.15 0.10 0.05 0.00

b 6

c

4 2

a

0

0.10 0.08 0.06 0.04 0.02 0.00

E10.5

1.6 1.2

ND5

a a

0.8 0.4 0.0

3.0

**

2.0 1.0 0.0 ND5

ND6 D-loop

Figure 6

IgG

5mC

**

15.0 **

10.0 5.0

**

**

**

**

0.0 **

60.0

E10.5

40.0 ** 20.0

**

**

IgG DNMT1

*

4.0

E7.5 20.0

0.0 12S

a

E3.5 E7.5 E10.5

D

E7.5

0.10 0.08 0.06 0.04 0.02 0.00

E6.5

2.0

E3.5 E7.5 E10.5

E7.5 E10.5

E

mt-Dnmt1

Methylation level detected by MeDIP-qPCR

Methylation level detected by MeDIP-seq

8

Relative expression levels

0.20

E3.5

C

Total Dnmt1

Relative expression levels

0.25

Fold enrichment

B RPKM of Dnmt1

A

**

**

ND6

A D

Nucleus MtGs for fusion and fission

Nrfs Polg, Polg2, Tfam, Ssbp1, Tk2, Stoml2, Fis1

Fis1, Mfn1, Mfn2, Opa1

Fis1

Opa1 Mfn1/2

E3.5

Fission

E7.5

Fusion

Fusion Fission Mitochondrial biogenesis & mtDNA replication

Casp3

B

Apoptosis

Damage

Cyto C

β-oxidation ROS TCA

Mitochondrion Slc25a11

Cyto C Bax, Bcl2l1, Bid, Bbc3 (Puma)

Damage

GSH

P53 OXPHOS

ROS

Nucleus ROS

GSH Nrf2 Kep1 Cul3

ROS

O-

O2-

Nrf2 Kep1 Cul3

Nrf2 Gs

Reduced glutathione (GSH) Gpx

Gclc, Gs, Gpxs, Gstm1, Glrx2

Gclc Glutamlycysteine

Clycine

Glutamate Cysteine

Gstm1 Oxidized glutathione (GSSG) ROS Prot-SH

Prot-S-SG Glrx2

Gclc

Upregulated genes

Mfn1

Downregulated genes Biochemical reaction Macromolecular damage Transmembrane transport Enhenced aerobic metabolism mtDNA

Figure 7

A

Ten most abundant MtGs in each cluster, enriched GO processes and KEGG pathways 10 most abundant MtGs

GO

Pathway

mt-Nd1, Hspa8, Hsp90ab1, Atp5l, mt-Cytb, mtNd2, mt-Co1, mt-Rnr2, mt-Rnr1, Cox5b

Oxoacid metabolism, Transmembrane transport, Translation, Electron transport chain, Coenzyme metaboliam

Oxidative phosphorylation, Huntington's disease, Alzheimer's disease, Parkinson's disease, Aminoacyl-tRNA biosynthesis

Cluster 2

Slc25a5, mt-Nd4, mt-Nd5, Ckb, Cox6a1, Ppp1cc, Cox8a, Uqcrc1, Calm1, Cs

Oxoacid metabolism; Fatty acid metabolism; Cellular lipid metabolism;Coenzyme metaboliam; Translation

Fatty acid metabolism, Alzheimer's disease, Propanoate metabolism, Valine, leucine and isoleucine degradation, Parkinson's disease

Cluster 3

Atp5b; Atp5g3; Atp5a1; Hspd1; Prelid1; Slc2a3; Vdac1; Atp5f1; Atp5e; Phb2

Oxoacid metabolism; Cellular amino acid metabolism, Transmembrane transport,Oxidative phosphorylation, Nitrogen compound biosynthesis

Parkinson's disease, Huntington's disease, Oxidative phosphorylation, Alzheimer's disease, Arginine, proline metabolism

Cluster 4

Rpl10a, Ywhae, Hspe1, Rpl23, Atp5h, Atp5g2, Ddx5, Cd24a, Uqcrh

Oxoacid metabolism; Apoptotic mitochondrial changes, Cellular amine metabolism; Mitochondrial transport; Cellular amino acid metabolism

Huntington's disease; Alzheimer's disease; Parkinson's disease; Arginine and proline metabolism; Oxidative phosphorylation

Cluster 5

Slc25a3, Ap2m1, Slc25a4, Atp5o Ndufb8, Atp5j, Ndufs3 Ndufb6, Gpx1, Mtx2

Electron transport chain; Response to reactive oxygen species; Transmembrane transport; Regulation of programmed cell death; Energy derivation by oxidation of organic compounds

Huntington's disease; Parkinson's disease; Alzheimer's disease, Oxidative phosphorylation; Insulin signaling pathway

Cluster 6

Stmn1, Cox7c, Sod1, Chchd2, Atp5j2, Idh2 Ndufa4, Ndufs2, Mdh1 Mtch1

Electron transport chain, Oxoacid metabolis, Fatty acid metabolism, Coenzyme metabolism, Cellular lipid metabolic process

Oxidative phosphorylation; Alzheimer's disease; Huntington's disease; Parkinson's disease; Valine, leucine and isoleucine degradation

B


Cluster 1

200 150

100 50 0

Supplementary Fig. 1. Dynamic and functional profiling of MtGs during early development. (A) Ten most

abundant MtGs in each cluster of Figure 1A, as well as enriched functional processes based on Gene ontology (GO) and KEGG analyses. (B) 3D-view of gene expression patterns of MtGs in each cluster of Figure 1A.

Supplementary Fig. 1