Mitochondrial DNA Deficiency and Supplementation in Sus scrofa Oocytes Influence Transcriptome Profiles in Oocytes and Blastocysts

Mitochondrial DNA (mtDNA) deficiency correlates with poor oocyte quality and fertilisation failure. However, the supplementation of mtDNA deficient oocytes with extra copies of mtDNA improves fertilisation rates and embryo development. The molecular mechanisms associated with oocyte developmental incompetence, and the effects of mtDNA supplementation on embryo development are largely unknown. We investigated the association between the developmental competence of Sus scrofa oocytes, assessed with Brilliant Cresyl Blue, and transcriptome profiles. We also analysed the effects of mtDNA supplementation on the developmental transition from the oocyte to the blastocyst by longitudinal transcriptome analysis. mtDNA deficient oocytes revealed downregulation of genes associated with RNA metabolism and oxidative phosphorylation, including 56 small nucleolar RNA genes and 13 mtDNA protein coding genes. We also identified the downregulation of a large subset of genes for meiotic and mitotic cell cycle process, suggesting that developmental competence affects the completion of meiosis II and first embryonic cell division. The supplementation of oocytes with mtDNA in combination with fertilisation improves the maintenance of the expression of several key developmental genes and the patterns of parental allele-specific imprinting gene expression in blastocysts. These results suggest associations between mtDNA deficiency and meiotic cell cycle and the developmental effects of mtDNA supplementation on Sus scrofa blastocysts.


Introduction
The mitochondrial genome (mtDNA) is a double stranded circular DNA that is approximately 16.6 kb in size and encodes 13 of the subunits of the electron transfer chain, which generates the vast majority of cellular ATP through oxidative phosphorylation (OXPHOS) [1]. It also encodes 2 rRNAs and 22 tRNAs and has one major non-coding region, the D-Loop, that is the site of interaction for the nuclear-encoded transcription and replication factors that translocate to the mitochondrion to drive first transcription and then replication [2]. Furthermore, mtDNA is inherited from the population present in the oocyte at the time of fertilisation and is, thus, a maternally-only inherited genome [3].
In human and other mammalian species, including the pig, there are associations between mtDNA copy number, low levels of ATP content, and the ability of the oocyte to fertilise and give rise to embryos [4][5][6][7][8]. For example, oocytes with low mtDNA are more likely to fail to fertilise or arrest during preimplantation development. Indeed, mitochondrial deficiency can lead to dysfunction in a number of cell types [9,10]. In the oocyte, mtDNA deficiency appears to be one of the causes of female factor infertility [4,5,8]. To overcome this problem, several approaches have been proposed and adopted. For example, oocytes can be supplemented with either ooplasm [11,12] or purified populations of mitochondria from autologous or  Cresyl Blue (BCB) labelling of cumulus-oocyte-complexes to assess cytoplasmic maturity of oocytes. After labelling of oocytes, BCB+ and BCB− oocytes were washed and cultured separately. Denuded metaphase II oocytes were pooled and used for RNAseq. (B) Schematic representation of the mtDNA supplementation procedure, i.e., fertilisation by intracytoplasmic sperm injection (ICSI) with (MT+) or without (MT−) extra copies of mtDNA. RNAseq data from BCB+ oocytes and individual blastocysts generated with or without mtDNA supplementation were used to perform longitudinal differential gene expression (DEG) analysis by comparing the transition from oocyte to either MT− or MT+ blastocysts. Details of the procedure to undertake the analysis are found in the Materials and Methods.
This study aimed to address two key questions using the pig as our model. First, we sought to determine the association between mtDNA deficiency and the transcription profiles of metaphase II oocytes. We undertook this utilising the most up to date version of the Sus scrofa reference genome and its annotations [35], and we included all the Ensembl annotated genes [36] to analyse single cell oocyte data [31]. Second, we sought to determine if the transcriptomic profiles during the developmental transition from the oocyte to the blastocyst would be affected following mtDNA supplementation. This was undertaken by employing a longitudinal analysis of differentially expressed genes [37,38]. Our analysis uncovered functional pathways affected in BCB− oocytes and altered the developmental gene expression in blastocysts as a result of mtDNA supplementation.

mtDNA Deficient Oocytes Revealed Downregulation of OXPHOS Genes
mtDNA deficiency in Sus scrofa oocytes correlates significantly with BCB staining [13,14]. We previously analysed RNAseq data sets comprising BCB positive (BCB+) and negative (BCB−) oocytes (Methods Section 4.1 and Figure 1A) to address the effect of mtDNA deficiency on global oocyte transcriptome profiles [31] using the Sus scrofa reference genome Sscrofa10.2 (Accession No. GCF_000003025.5, released on 7 September 2011) [32]. However, more than half of the differentially expressed genes (DEGs) did not have proper annotation with gene symbols, and the results did not include read counts for the mtDNA encoded genes [31]. We reanalysed these BCB+ and BCB− oocyte data sets (Table S1) using the latest reference genome Sscrofa11.1 Accession No. GCF_000003025.6 [35] and Ensembl annotation release 105 [36] to determine if there would be greater clarity about the genes affected and the associated functional pathways and networks.
Single oocyte RNAseq data contained high numbers of sample replicates, i.e., 15 BCB+ and 14 BCB− oocytes (Methods Sections 4. 1 and 4.4). In all, 10 to 30 million reads per sample were analysed, and the levels of read mapping and assignment of annotation (Table  S1) were comparable with previous analyses [39,40]. Principal component analysis (PCA) revealed high consistency for oocyte RNAseq profiles relative to blastocyst RNAseq data, as previously shown [19] ( Figure S1A). We found that the RNAseq profiles of the oocytes were significantly affected by ovary source ( Figure S1B, left panel) suggesting a significant effect of genetic background on each oocyte's transcriptome profile. Therefore, ovary source was included in the data analysis, and the effect of batch was removed ( Figure S1B, right panel).
An analysis of the transcriptome profiles between BCB+ and BCB− oocytes identified 402 DEGs at a significance of FDR < 0.05 (Table 1 and Table S2). Amongst those, mtDNA encoded genes were abundant, including 11 protein coding genes and five tRNAs that were all significantly downregulated in BCB− oocytes ( Figure 2A). In addition, cytochrome c oxidase I (COXI) and NADH dehydrogenase 3 (ND3) also showed downregulation at an FDR < 0.1 indicating all protein coding genes of mtDNA were downregulated in BCB− oocytes. This is consistent with the lower levels of mtDNA copy number found in BCB− oocytes [13] resulting in reduced transcriptional activity for the mitochondrial genome. Gene ontology (GO) enrichment analysis of the whole DEG data sets identified genes involved in OXPHOS including ATP synthesis coupled electron transport, oxidative phosphorylation, and mitochondrial ATP synthesis, which were enriched (Table S3). Therefore, mtDNA deficiency in BCB− oocytes not only affected mtDNA encoded gene transcription, but also resulted in the downregulation of the majority of OXPHOS genes, which includes genes encoded by the nucleus ( Figure 2B and Figure S2).

RNA Metabolism and Meiosis/Mitosis Related Gene Expression Were Affected in BCB− Oocytes
DEG analysis also identified 308 protein coding nuclear genes as well as 74 noncoding transcripts (Table S2). Ten out of the top 30 DEGs were small nucleolar RNAs (snoRNA) ( Table 1) with a total of 56 snoRNA found to be downregulated in BCB− oocytes ( Table 2). The primary function of snoRNA is the chemical modification of RNAs, mainly ribosomal RNA, transfer RNA, and small RNAs [41,42]. A recent study identified a noncanonical function for snoRNA, such as the regulation of chromatin structure and mRNA splicing. Consequently, the downregulation of a large subset of snoRNAs could affect RNA metabolism. Indeed, genes associated with mRNA processing and metabolism were enriched in the DEGs (Table S2), and therefore, corresponding GO and REACTOME terms (Tables S4 and S5) were identified as functional network groups ( Figure 3B).    tology (GO) enrichment analysis of the whole DEG data sets identified genes involved in OXPHOS including ATP synthesis coupled electron transport, oxidative phosphorylation, and mitochondrial ATP synthesis, which were enriched (Table S3). Therefore, mtDNA deficiency in BCB− oocytes not only affected mtDNA encoded gene transcription, but also resulted in the downregulation of the majority of OXPHOS genes, which includes genes encoded by the nucleus (Figures 2B and S2). The entire OXPHOS related protein interaction network is shown in Figure S2. Levels of differential gene expression between BCB+ and BCB− oocytes are shown as fold change by colour scale. Red and blue represent up-and down-regulation in BCB− oocytes. Genes in grey box have no DEG data due to low or no expression.  splicing. Consequently, the downregulation of a large subset of snoRNAs could affect RNA metabolism. Indeed, genes associated with mRNA processing and metabolism were enriched in the DEGs (Table S2), and therefore, corresponding GO and REACTOME terms (Tables S4 and S5) were identified as functional network groups ( Figure 3B).  Table S4.   Table S4.
The most affected functional network in BCB− oocytes was the meiosis and mitosis cell cycle process ( Figure 3A). A total of ten groups of GO terms associated with oocyte maturation and processes of cell division were found to be enriched in DEGs (Table S4). This suggests that poor oocyte quality is associated with the mis-regulation of genes required for the completion of metaphase II and following the first zygotic cell division after fertilization with the male gamete. Overall, mtDNA deficiency in BCB− oocytes affects the transcription of genes encoded by the nuclear genome that are essential for oocyte development and maturation and RNA metabolism, which, in turn, could result in the oocyte being defective.

The Use of mtDNA Supplementation in Oocytes Influences the Expression of Genes Involved in Development and Differentiation in Blastocysts
Transcriptional profiles change greatly during the developmental transition from a quiescent state in the oocyte to an active state in the developing embryo [43,44]. We have previously demonstrated the use of longitudinal analysis by assessing differentially methylated regions of the Sus scrofa genome from oocytes to blastocysts, which provided an additional layer of information and insights into the effect of mitochondrial supplementation [19]. We utilized this approach to analyse our RNAseq data sets of oocytes and blastocysts generated by intracytoplasmic sperm injection accompanied by (MT+) or without mtDNA (MT−) supplementation ( Figure 1B). We hypothesized that the longitudinal analysis [37,38] would highlight the effect of mtDNA supplementation on blastocyst generation. To do this, DEGs between oocytes and blastocysts without mtDNA supplementation and between oocytes and blastocysts with supplementation were identified independently (Methods Sections 4.1-4.4 and Figure 1B). We then compared the two lists of DEGs to determine if there were any unique DEGs in either of the developmental transition processes ( Figures 1B and 4A).
Approximately 11,000 genes were differentially expressed between the oocyte stage and blastocyst stage, of which 90% of DEGs (10,000 genes) were common irrespective of whether mtDNA supplementation had taken place ( Figure 4A). In all, 891 genes were identified from the list of DEGs for the transition from oocyte to blastocyst generated without additional mtDNA. Of those genes, 136 had an FDR < 0.01 (Table 3 and Table S6). However, this subset of genes was not significantly different between the transition from oocytes to blastocysts that were supplemented ( Figure S3A). Similarly, 1547 genes were only identified in the DEGs from the transition from oocyte to blastocysts supplemented with mtDNA. Of these, 475 genes had an FDR < 0.01, and most were highly expressed in the blastocysts ( Figure S3B and Table S7). Some of these genes have been previously identified with a direct comparison between supplemented and non-supplemented blastocysts through RNAseq. For example, these include several tRNAs, sestrin 1 (SESN1), glycerophosphocholine phosphodiesterase 1 (GPCPD1), and NIMA related kinase 2 (NEK2) [19]. Many others were successfully uncovered through the longitudinal analysis, complementing the direct comparison method, and further addressing the effect of mtDNA supplementation in a developmental context.
The non-supplemented blastocysts showed slightly less transcriptional change during transition and the relatively lower number of DEGs were unique to this transition ( Figure 4A). Only four associated GO and REACTOME terms were enriched (Table S8). On the other hand, mtDNA supplementation of oocytes induced greater levels of gene activation in blastocysts, mainly nuclear-encoded protein coding genes (Table S7) and 211 GO and REACTOME terms, which consisted of 53 functional annotation groups that were enriched in this subset of DEGs (Table S9). Genes involved in cell differentiation and various tissue development, including the pancreas, glands, and neurons, were abundant and uniquely regulated in the process of transition to blastocyst development following supplementation ( Figure 4B,C). Several of these genes, such as BCL2 antagonist/killer 1 (BAK1), neuregulin 1 (NRG1), platelet derived growth factor receptor alpha (PDGFRA), and hepatocyte nuclear factor 4 alpha (HNF4A), have critical roles in developmental processes [45][46][47][48]. The lack of activation of these genes in the non-supplemented blastocysts might have detrimental effects on early embryo development or be indictive of the premature expression of genes of early differentiation. Overall, the results of the longitudinal DEG analysis highlighted a novel aspect of the oocyte-to-blastocyst transition process.  Network associated with gland and pancreas development. GO biological process and REACTOME terms in the same group are coloured and presented with their representative group name. The number for each group name corresponds to the GO and REACTOME group number in Table S9.  (C) Network associated with gland and pancreas development. GO biological process and REAC-TOME terms in the same group are coloured and presented with their representative group name. The number for each group name corresponds to the GO and REACTOME group number in Table S9.

Effect of mtDNA Supplementation on Parent of Origin Gene Expression Patterns
We have previously shown that mtDNA supplementation was associated with only minor differences in the DNA methylation of imprinted genes in Sus scrofa, and no differences were found in the levels of imprinted gene expression between supplemented and non-supplemented blastocysts [19]. However, it is still unknown if they showed parent of origin imprinted gene expression patterns. Therefore, we investigated the allele-specific expression of imprinted genes and neighbouring genes by analysing the frequency of single nucleotide polymorphisms (SNP) [49]. Since we used commercial pig sperm and oocytes collected from ovaries processed in a local Australian abattoir, it was not feasible to track parental genotypes from available resources. Therefore, we analysed SNPs in the transcripts, which reflect the expression of both parental alleles. If only one of the parental alleles is dominantly active for transcription as in imprinted genes, we would expect to find no or low numbers of SNPs in the transcript.
The method requires genetic variations between the two parental pigs. Average SNP frequency in the human genome is about one in every 300 base pairs with a minor allele frequency greater than 1% [50][51][52], and an NGS analysis of four swine breeds using reduced representation genomic libraries identified 372K SNPs [53]. Therefore, it seems feasible that the allele-specific expression could be analysed by SNP identification with the depth of our blastocyst RNAseq data. We set up detection criteria for the number of reads per transcript, depth of sequence reads, number of SNPs per transcript, and level of a variant, as described in the Methods, to analyse allele-specific expression of genes at the imprinting loci. For example, insulin like growth factor 2 receptor (IGF2R) on Sus scrofa chromosome 1, one of the most characterized imprinted genes in human and mouse genomes [54][55][56], encodes 7560 bp transcripts (Accession no NM_001244473). In blastocyst sample MT+ #5, derived through supplementation, its RNAseq data identified three SNPs with a mean variant frequency of 52% ( Figure 5A). In the same data, 27 SNPs were identified in a total transcript length of 29,281 bp containing 10 neighbouring genes. SNP frequency, calculated by SNP number divided by transcript length (kb), was used as an indicator of mono-allelic or bi-allelic expression at low and high frequency, respectively. In the RNAseq data set MT+ #5 (mtDNA supplemented blastocyst), the SNP frequency for IGF2R was 0.40 SNP per kb transcript (SNP/kb) whilst there was a 0.92 SNP/kb frequency for its neighbouring genes. This suggests that IGF2R was biased to mono-allelic expression as expected for imprinted genes. Similarly, in RNAseq data set MT− #3 (non-supplemented), the SNP frequency for IGF2R was 1.98 SNP/kb and 0.61 SNP/kb for its neighbouring genes. The even higher SNP frequency for IGF2R than its neighbouring genes suggests that both parental alleles were actively transcribed for this gene in this blastocyst.
It is possible that a low or high SNP frequency may simply reflect the genetic variation of parents for specific genes and genomic regions. To minimize distribution bias of the SNPs in specific genomic regions, which affects bi-allelic expression, we analysed 10 imprinted genes and 73 non-imprinted neighbouring genes at six imprinting loci across four chromosomes. Total transcript lengths of 57,219 bp and 273,582 bp were analysed for imprinted genes and non-imprinted neighbouring genes, respectively, and we normalized SNP data by length for comparison. Table S10 summarises the SNP analysis results, and the sex of each blastocyst was included in the model of statistical analysis. For genes with bi-allelic expression, about 50% variant frequency is expected, and we found~40% variant frequency in non-imprinted genes from both supplemented and non-supplemented blastocysts ( Figure 5B). Variant frequency in imprinted genes was similar, and this was understandable as we filtered and kept SNPs with the signature of bi-allelic expression. As mentioned above, it would be SNP frequency (SNP number per transcript length) that provides the status of allelic expression. In MT− blastocysts, there were no differences in SNP frequency between imprinted and non-imprinted genes ( Figure 5C, left panel). This was interesting as the mono-allelic expression of imprinted genes is expected to contribute to a lower SNP frequency more than non-imprinted genes. Indeed, supplemented blastocysts showed a significantly lower SNP frequency in imprinted genes than non-imprinted genes ( Figure 5C, right panel). As we analysed over 330 kb total transcript length for 83 genes, these were unlikely to be caused by specific genotype combinations of specific genes in non-supplemented blastocysts. However, this might suggest that blastocysts generated without mtDNA supplementation have delayed or disturbed imprinted gene expression or partial inactivation of one of the parental genomes for regular bi-allelic expression, or a combination of both. It is possible that a low or high SNP frequency may simply reflect the genetic var tion of parents for specific genes and genomic regions. To minimize distribution bias the SNPs in specific genomic regions, which affects bi-allelic expression, we analysed imprinted genes and 73 non-imprinted neighbouring genes at six imprinting loci acro four chromosomes. Total transcript lengths of 57,219 bp and 273,582 bp were analysed imprinted genes and non-imprinted neighbouring genes, respectively, and we norm ized SNP data by length for comparison. Table S10 summarises the SNP analysis resu and the sex of each blastocyst was included in the model of statistical analysis. For gen

Discussion
In this study, we have analysed RNAseq data sets from Sus scrofa oocytes and blastocysts generated by injecting sperm into a mature oocyte with or without mtDNA supplementation. RNAseq data were analysed globally using tagwise dispersion estimation and Generalized Linear Models (Methods Section 4.4) which provide unbiased estimates of the gene-specific dispersions. This approach is suitable for replicates with high variation, such as blastocyst RNAseq data. We thought validation by other tests, such as qRT-PCR, would be difficult for ultra-low quantity materials considering the associated technical and handling errors, selection of appropriate controls and detection thresholds, and variations amongst samples. The use of the most up-to-date Sus scrofa reference genome sequence and associated annotations identified differentially expressed genes in BCB− oocytes, including a large subset of snoRNAs as well as genes associated with meiotic and mitotic cell cycle processes, that were not identified in our previous analysis [31]. Furthermore, using a longitudinal transcriptome analysis approach [37,38] likely provided a greater range of DEGs and revealed activation of subsets of developmental genes in mtDNA supplemented blastocysts. This approach circumvented the issue of fewer genes being identified as a result of the direct comparison between MT+ and MT− blastocysts due to the high FDR values associated with high sample variations [19].

mtDNA Deficiency in BCB− Oocytes Associates with Downregulation of OXPHOS Genes
The developmental competence of oocytes has been assessed with BCB labelling in various mammals with BCB− oocytes reportedly associated with a lower fertilisation potential, lower frequency of development to blastocyst, lower levels of mitochondrial activity, changes in expression of genes involved in apoptosis, and lower levels of mtDNA copy number [21,24,25,28,30,57]. Our results provide further evidence that BCB− oocytes have highly distinct transcription profiles compared to BCB+ oocytes (Table 1 and Table  S2 and Figure 2). A total of 16 mtDNA encoded genes, including 11 protein coding and 5 t-RNAs, were significantly downregulated in BCB− oocytes with two mtDNA protein coding genes, COXI and ND3, also downregulated but to levels of lower significance (FDR < 0.1). Since mtDNA copy number in pig BCB− metaphase II oocytes is~3-fold lower than in BCB+ oocytes [13] and that there are strong associations in other species, such as the human [8], pig [58], and mouse [59], with mtDNA copy number and fertilisation outcome, this was not particularly surprising. However, it confirmed low levels of expression for all of the mtDNA-encoded genes that are associated with OXPHOS and that the majority of nuclear-encoded genes associated with OXPHOS were also downregulated in BCB− oocytes (Table 1 and Table S2 and Figure 2 and Figure S2). On the other hand, we did not find differences in the expression of genes involved in mtDNA replication, such as DNA polymerase gamma (POLG), mitochondrial transcription factor A (TFAM), and nuclear respiratory factor 1 (NRF1), which have been previously reported to be differentially expressed in BCB− oocytes [30,57], although it should be pointed out that the previous studies analysed these genes in isolation with qRT-PCR and not as part of large data sets. Overall, the association between BCB− oocytes and the downregulation of OXPHOS-related genes are consistent with the previous observations that BCB− oocytes have lower levels of mtDNA copy number and mitochondrial activity [13,21,24,25].

Downregulation of RNA Metabolism and Meiotic and Mitotic Cell Cycle Process Genes in mtDNA Deficient Oocytes
One of the most striking differences between the transcriptome profiles for BCB+ and BCB− oocytes was the downregulation of a large subset of snoRNAs (Table 2). snoRNAs play a central role in ribosome biogenesis, guiding the sequence-specific chemical modification of pre-rRNA and its processing [41,42,60]. snoRNAs are components of ribonucleoprotein complexes and are divided into two main classes based on highly conserved sequences, namely box C/D and box H/ACA motifs. In addition, some snoRNAs of both classes are localized in Cajal bodies, known as small Cajal body-associated RNAs (scaRNA). We determined these three types of snoRNAs were downregulated in BCB− oocytes ( Table 2). The main functions of canonical snoRNAs and scaRNAs are their catalysation of 2'-O-ribose methylation (box C/D) and guidance of the pseudouridylation (box H/ACA) of rRNA and small nuclear RNAs, respectively. Subsets of snoRNAs have noncanonical functions, including the regulation of splicing, mediators of oxidative stress, and regulation of chromatin structure [41,61]. In vertebrates, snoRNAs are encoded in the introns of host genes, protein coding mRNA or long noncoding RNA, and released from the host gene precursor RNAs by splicing. After the splicing and production of snoRNAs, host gene RNAs are directed to degradation via nonsense-mediated RNA decay (NMD) [62], and NMD genes were significantly affected in BCB− oocytes (REACTOME Group 3 in Table S5). Although there are known functions for snoRNAs in RNA biogenesis and metabolism, as mentioned above, the functional annotations of snoRNA genes (SNORD and SNORA) themselves were not readily available for the GO and REACTOME enrichment analysis [63,64]. Therefore, snoRNAs did not contribute to the enrichment analysis results, and no SNORD and SNORA genes were found in Tables S3-S5 (in 'Associated Genes' columns). Despite this, the RNA metabolic process, catabolic process, and splicing and processing were still found to be significantly enriched (Tables S4 and S5 and Figure 3B) indicating that RNA metabolism is largely affected in BCB− oocytes.
The most affected functional network in BCB− oocytes was meiotic and mitotic cell cycle processes, including chromosome segregation, spindle organization, DNA replication, centrosome duplication, and cell cycle phase transition ( Figure 3A and Tables S4 and S5). After maturation in the ovary or in medium, oocyte development is arrested at the metaphase II stage prior to fertilisation [65]. Fertilisation triggers the completion of metaphase II by segregating half of the sister chromatids into the second polar body and forming a female haploid pronucleus to fuse with the male pronucleus. The zygote then assembles the first mitotic spindle for the first cell division. Genes involved in these cell cycle processes were differentially expressed in BCB− oocytes and were mostly downregulated (Tables S2 and S5). It is reasonable to assume that the lack of these critical gene transcripts in readiness for the completion of meiosis II and following the first mitotic division could cause significant developmental defects at fertilisation and during embryogenesis. These might include failed chromosomal segregation and subsequent aneuploidies that lead to embryonic arrest or spontaneous abortion [66,67]. The downregulation of RNA metabolism and processing is likely to be part of the cause of the under-representation of those critical cell cycle gene transcripts.

mtDNA Supplementation Maintains Expression of Genes Necessary for Embryogenesis
Longitudinal DEG analysis identified subsets of genes that are differentially expressed in mtDNA supplemented-derived blastocysts, developed from quiescent state oocytes (Table S7). Those genes were mostly upregulated in the supplemented blastocysts and involved in various developmental processes (Figure 4 and Table S9). These include BAK1 which functions as a pro-apoptotic regulator involved in a wide variety of cellular activities and shaping tissue and organ morphologies via apoptosis [48,68,69]. In addition, PDGFRA is a receptor kinase binding to PDGF that stimulates the growth and migration of vascular smooth muscle cells, fibroblasts, and glial cells [46,70,71]. On the other hand, deletion of PDGFRA induces embryonic lethality and causes the defective development of many endoderm-and mesoderm-derived structures [72,73]. Zinc finger protein 792 (ZNF792) is a paralogue of ZNF304, which recruits repressive cofactors and induces DNA hypermethylation and transcriptional silencing [74]. HNF4A is a nuclear transcription factor that regulates the expression of several hepatic genes for the development of the liver, kidney, and intestines. Mutations in HNF4A are associated with a monogenic form of type 2 diabetes and the loss of function in mice, causing severe hepatomegaly and steatosis and resulting in premature death [45,75]. Moreover, NRG1 is a trophic factor that contains an epidermal growth factor (EGF)-like domain and NRG1 signalling is required for embryonic development. The loss of function of NRG1 can cause homozygous lethality and behavioural abnormalities in heterozygous mutant mice [47,76] whilst the addition of NRG1 to the oocyte in vitro maturation media results in improved embryo development [58]. These genes were activated in the supplemented blastocyst during the oocyte-to-blastocyst transition process, but not significantly upregulated in the non-supplemented blastocysts. Since mtDNA supplementation induces mtDNA replication in early embryogenesis and improves embryo quality [13], this likely enhances developmental gene expression during progression to the blastocyst stage.
Communication between the nuclear and mitochondrial genomes is important for effective cellular function. It is well-documented that ATP generated through OXPHOS requires the replication of the mitochondrial genome to have taken place [77] and that mtDNA copy is not uniform across all cell types but rather reflects the requirements of each cell type for OXPHOS-derived ATP [78]. Indeed, throughout development, there are changes in mtDNA copy that reflect the metabolic requirements of the developing embryo and foetus. For example, undifferentiated (pluripotent) cells have a requirement for low levels of OXPHOS-derived ATP [79], as is the case for tumour-initiating cells [80] whilst heart cells and neurons have high levels of the mtDNA copy number to support their complex functions [81]. A failure to mediate the interactions between the two genomes will inhibit key developmental processes, such as differentiation, and result in cellular arrest. However, the modulation of one genome can impact the other. For example, the use of DNA demethylation agents on tumour-initiating cells can result in increased differentiation potential and increased mtDNA copy number that take place in a synchronous fashion [82]. Likewise, the short-term depletion of mtDNA copy number can reset the differentiation potential of tumour-initiating cells whilst longer term depletion can inhibit tumorigenesis with both outcomes exhibiting differences in nuclear gene expression [80]. mtDNA deficiency in oocytes likely provides another example of how one genome is asynchronous with the other. Indeed, BCB− oocytes are indicative of poor developmental competence where nuclear and cytoplasmic maturation are not in synchrony. As a result, these oocytes frequently either fail to fertilise or arrest during preimplantation development [14,30], but they can be rescued through mtDNA supplementation [13]. Developmental competence can also be rescued through the addition of agents like NRG1 [58] or the glycoside mogroside V [83] to the oocyte maturation media, which impact the mtDNA copy number and enhance embryo development.

mtDNA Supplementation Might Improve Imprinted Gene Expression
Accumulated evidence suggests that assisted reproductive technologies, for example, in vitro embryo culture, in vitro fertilisation, and intracytoplasmic sperm injection affect imprinted gene methylation and expression compared to the children from spontaneous conception [84][85][86][87]. Embryos generated by in vitro fertilisation and cultured in medium can either show the bi-allelic expression of H19 or limited paternal expression [85,88]. Moreover, differential DNA methylation and expression of IGF2/H19 and IGF2R in placenta and cord blood have been reported for in vitro conceived children [89]. Furthermore, altered expression of H19 and Pleckstrin homology-like domain family A member 2 (PHLDA2) has been observed in human placentas following in vitro fertilisation and intracytoplasmic sperm injection, but this was not due to the loss of the pattern of imprinted gene expression [90]. Our previous study identified differentially methylated regions throughout the nuclear genome between blastocysts derived through either mtDNA supplementation or without, but they showed only minor differences in DNA methylation in imprinting control regions (ICRs) [19]. Differential DNA methylation patterns of ICRs, inherited from paternal and maternal genomes, are essential for imprinted gene expression and maintained in somatic cell lineages [55]. Although there was no difference in the level of imprinted gene expression associated with mtDNA supplementation [19], we observed a potential delay or disturbance of the establishment of the patterns of parental allele-specific imprinted gene expression in the non-supplemented blastocysts ( Figure 5 and Table S10). We thought the effect of blastocyst sample variation on the imprinted gene expression pattern was minimal as DNA methylation of ICRs was already established prior to fertilisation and maintained in the blastocysts [19]. Blastocysts derived from mtDNA supplementation showed less bi-allelic expression in imprinted genes and higher levels in neighbouring non-imprinted genes, suggesting an improved establishment of parental allele-specific expression patterns. Combining mtDNA supplementation with specific embryo culture conditions [91] might be a potential approach to improve gene expression profiles so that they more closely resemble the profiles of naturally conceived embryos. However, this would require a full understanding of how the mtDNA supplementation approach enhances embryo development.
We have shown that mtDNA deficiency in Sus scrofa oocytes significantly affects the expression of genes involved in RNA metabolism and processing. Furthermore, the downregulation of genes associated with meiotic and mitotic cell cycle processes are likely to be one of the direct causes of poor oocyte quality and poor fertilisation outcome. However, further understanding of how mtDNA deficiency influences nuclear gene expression with respect to communication between the two genomes needs to be established. We also showed that the mtDNA supplementation of oocytes in combination with sperm injection improves the maintenance of developmental gene expression and patterns of parental allele-specific imprinted gene expression in blastocysts. These results further support the concept that the addition of extra copies of mtDNA has a significant impact on the overall gene expression profiles during oocyte-to-blastocyst transition and contribute to improved fertilisation and preimplantation development outcomes.

Oocyte Collection and BCB Staining
Cumulus-oocyte complexes (COCs) were collected from gilt ovaries obtained from a local abattoir and cultured as previously described [19,31]. Briefly, porcine ovaries were transported to the laboratory in warm 0.9% NaCl solution. The COCs were then aspirated from follicles with diameters of between 3 and 6 mm using an 18 G needle. The COCs were then washed 3 times in handling media (25 mM Hepes-TCM199, Gibco ® , Waltham, MA, USA) supplemented with 10% sow follicular fluid (SFF). The COCs were incubated in 12 µM BCB in in vitro maturation (IVM) medium for selection as described in [13], and BCB+ and BCB− COCs were then washed and cultured separately for 42-44 h in pre-equilibrated IVM medium. To collect metaphase II (MII) oocytes, expanded COCs following IVM were briefly treated with 0.1% (0.5mg/mL) of hyaluronidase and gentle pipetting. The resultant denuded oocytes were washed using a narrow glass pipette to completely remove all cumulus cells ( Figure 1A). Cohorts of single MII oocytes were collected and stored at −80 • C prior to RNA extraction [31].

Generation of Blastocysts by Intracytoplasmic Sperm Injection (ICSI) with (MT+) or without (MT−) mtDNA Supplementation
Expanded blastocyst stage Sus scrofa embryos were generated from in vitro matured oocytes by ICSI; and mitochondrial supplementation in combination with ICSI using third-party oocytes (heterologous) or sister oocytes (autologous) as the source of the mitochondrial isolate as previously described in [13,19]. Briefly, 15 oocytes were homogenised in mitochondrial isolation buffer (20 mM Hepes pH 7.6, 220 mM mannitol, 70 mM sucrose, 1 mM EDTA) containing 2 mg/mL BSA by a drill-fitted Teflon pestle. The homogenate was centrifuged at 800× g for 10 min at 4 • C to remove cell debris, and then, supernatants were centrifuged at 10,000× g for 20 min at 4 • C to pellet the mitochondrial fraction. Following two washing steps, the mitochondrial pellet was resuspended in 5 µL of mitochondrial buffer. For ICSI with mitochondrial supplementation, 3 pl of mitochondrial isolate along with a single sperm was injected into the oocyte as described [13]. Single embryos were stored in 0.2 mL PCR tubes in a −80 • C freezer prior to RNA extraction.

RNA Extraction, RNAseq Library Preparation, and Sequencing
Total RNA was extracted from single oocytes or single expanded blastocysts using the PicoPure RNA isolation kit (Thermo Fisher Scientific, Waltham, MA, USA), according to the manufacturer's instructions. The quality of total RNA was assessed using High Sensitivity RNA ScreenTape (Agilent Technologies, Santa Clara, CA, USA). Next generation sequencing (NGS) libraries were prepared using Ovation ® RNA-Seq system V2 and Ovation ® Ultralow System V2 (NuGEN Technologies, Männedorf, Switzerland) for single MII oocytes [31] and the Trio RNA-Seq Library Preparation Kit (Tecan Group Ltd., Männedorf, Switzerland) for single blastocysts [19]. NGS libraries were sequenced using 100 bp paired-end sequencing chemistry.

RNAseq Data Analysis and Differentially Expressed Gene (DEG) Identification
Single oocyte (NCBI Sequence Read Archive IDs: SRR6451030-SRR6451058) and blastocyst (SRR16706437-SRR16706452) RNAseq data from our previous studies [19,31] were used. Oocyte (14 BCB− oocytes and 15 BCB+ oocytes) and blastocyst (5 MT− and 5 MT+) RNAseq sample metadata can be found in Table S1. RNAseq raw fastq files were quality checked with 'fastqc' (version 0.11.9) [92], and trimming of adapters and quality filtering were then performed with 'fastp' (version 0.20.1) [93] with options: -detect_adapter_for_pe, -q 20, -length_required 30. Trimmed and quality filtered paired-end reads were aligned to the Sus scrofa reference genome sequence Sscrofa11.1 accession No. GCF_000003025.6 [35] and Ensembl annotation release 105 [36] using 'STAR' (version 2.7) [94] with default parameters. Gene expression was quantified by counting the number of reads aligned to each Ensembl gene model using 'featureCounts' (version 1.5.2) [95], and output results were assessed for mapping quality by MultiQC version 1.9 [96]. Summary statistics for the RNAseq data are shown in Table S1. The Trimmed Mean of M-values (TMM) normalisation method from edgeR was applied to normalise read counts according to library size differences between samples. PCA was performed to visualise the summary of gene expression profiles for all RNAseq data. We found that the source of the ovary affected gene expression profiles identified by PCA; therefore, this was included in the linear mixed model as a covariate and batch effect was successfully corrected in the expression data.
For comparisons between BCB+ and BCB− oocytes, genes with low expression were filtered out prior to DEG analysis, keeping genes with at least 1 count per million (CPM) reads in more than eight samples. Data were analysed by using generalized linear models 'glmfit' and 'glmLRT' in edgeR version 3.32.0 [97] in R version 4.0.3 with the source of ovary as a covariate factor in the model. Genes were considered differentially expressed if their FDR (false discovery rate) or adjusted p-value was <0.05 (Table S2). Similarly for longitudinal DEG analysis, DEG analysis was carried out on RNAseq data from BCB+ oocytes and blastocysts generated with (MT+) or without (MT−) mtDNA supplementation, independently ( Figure 1B). Then, two subsets of DEGs were compared to identify unique DEGs for the transition from oocyte to either MT− or MT+ blastocysts. In this analysis, genes with low expression were filtered out prior to DEG analysis, keeping genes with at least 1 count per million (CPM) reads in more than four samples. DEGs between BCB+ oocyte and MT− or MT+ blastocysts were identified using the limma-voom method (version 3.46.0) [98,99] to maintain consistency with the analysis tool used in a previous report [19]. An FDR threshold < 0.1 was used for significant DEGs between BCB+ oocytes and blastocysts with or without mtDNA supplementation and compared two subsets of DEGs for the efficient elimination of common DEGs. A Venn diagram was generated to show common and unique subsets of DEGs for each comparison ( Figure 4A), and high confidence (FDR < 0.01) unique DEGs are listed in Tables S6 and S7. Results were visualized with jittered-boxplot, using ggplot2 [100].

Functional Pathway Enrichment and Gene Network Analysis
Functional pathway enrichment analyses for DEGs were performed using the GSEA v4.2.3 [64] and Enrichment Map [101] as described in [102]. Briefly, Sus scrofa Ensembl gene IDs were used to search for corresponding human orthologue gene symbols using the Ensembl BioMart database [103] and used for pathway enrichment analysis input. Annotation gene set file, Human_GO_AllPathways_with_GO_iea_May_25_2022_symbol.gmt from Bader lab gene sets collections (http://download.baderlab.org/EM_Genesets/current_ release/, accessed on 8 February 2023), was used for GSEA analysis with a default FDR threshold of 0.25 (Table S3). The enrichment analysis results for GO biological process [104] and REACTOME pathway [105] were then visualized with the Enrichment Map app in Cytoscape [106]. We also used ClueGO [63] for pathway enrichment analysis with subsets of DEGs (Tables S4, S5, S8 and S9). For protein-protein interaction network analysis, the STRING [107] app in Cytoscape was used.

Bi-Allelic Expression Analysis of Genes in Imprinting Loci
Imprinted genes and neighbouring genes at the locus were analysed for the level of bi-allelic expression by identifying single nucleotide polymorphisms (SNP) in the mtDNA supplemented (MT+) and non-supplemented (MT−) blastocyst transcripts. Sus scrofa imprinted genes (https://www.geneimprint.com/site/genes-by-species.Sus+scrofa, accessed on 8 February 2023) with sufficient read depth were selected for analysis. We set a total of 5000 reads or more per imprinted gene transcript from ten blastocyst RNAseq data sets (MT−: n = 5 and MT+: n = 5) as a threshold to obtain sufficient read counts for SNP identification (Table S10). We analysed 10 imprinted genes and 73 non-imprinted neighbouring genes at six imprinting loci across four chromosomes. Imprinted genes and loci selected for investigation were: Chr1: IGF2R and NDN loci, Chr6: OSBPL1A and DIRAS3 loci (including SGIP1), Chr9: PEG10 locus (including CASD1, PPP1R9A), Chr14: INPP5F locus (including TACC2).
SNPs were identified and counted using 'samtools mpileup' for mapped reads pileup [108] and 'VarScan' version 2.3.8 [109] for variant sequence search. Imprinted gene locus sequences, including 2Mb upstream and downstream sequences of the Sus scrofa reference genome sequence Sscrofa11.1 accession No. GCF_000003025.6, were used for read pileup. Variant call by VarScan was made with the options: -min-avg-qual 20 -min-var-freq 0.01 -min-reads2 10, and the results were further filtered by keeping variant calls if there were >25 read counts per site and variant frequency was between 20-80%. Furthermore, if there were more than three SNPs per transcript, they were kept for further analysis as confident SNPs indicating bi-allelic expression status. SNPs less than three per transcript were discarded as they might be associated with technical issues, such as PCR and sequencing errors. SNP frequency was calculated by number of SNP counts per kb of transcript and mean variant frequency was also calculated (Table S10). ANOVA test was carried out to see statistical significance associated with blastocyst type (MT− and MT+) and imprinting status by including sex of blastocyst, date of sampling, and date of RNAseq library construction as covariate factors. Sex of each blastocyst was determined by assessing read counts of selected chromosome X-and Y-linked genes [110].

Data Availability Statement:
The datasets supporting the conclusions of this article are available in the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra, accessed on 8 February 2023) under the two BioProject IDs PRJNA429045 and PRJNA777282.