Understanding genetic incompatibilities and genetic introgression between incipient species are major goals in evolutionary biology. Mitochondrial genes evolve rapidly and exist in dense gene networks with coevolved nuclear genes, suggesting that mitochondrial respiration may be particularly susceptible to disruption in hybrid organisms. Mitonuclear interactions have been demonstrated to contribute to hybrid dysfunction between deeply divergent taxa crossed in the laboratory, but there are few empirical examples of mitonuclear interactions between younger lineages that naturally hybridise. Here we use controlled hybrid crosses and high resolution respirometry to provide the first experimental evidence in a bird that inter-lineage mitonuclear interactions impact mitochondrial aerobic metabolism. Specifically, respiration capacity of the two mitodiscordant backcrosses (with mismatched mito-nuclear combinations) differ from one another, although they do not differ significantly from the parental groups or mitoconcordant backcrosses as we would expect of mitonuclear disruptions. In the wild hybrid zone between these subspecies the mitochondrial cline centre is shifted west of the nuclear cline centre, which is consistent with the direction of our experimental results. Our results therefore demonstrate asymmetric mitonuclear interactions that impact the capacity of cellular mitochondrial respiration and may help to explain the geographic discordance between mitochondrial and nuclear genomes observed in the wild.

Study species

Long-tailed finches of each subspecies (Poephila acuticauda acuticauda and P. a. hecki) used in this study came from captive populations maintained at Macquarie University in Sydney, Australia. Wild caught P. a. acuticauda individuals were originally sourced from Mount House (17°02′S, 125°35′E) and Nelson’s Hole (15°49′S, 127°30′E) in Western Australia, and wild-caught P. a. hecki were sourced from October Creek (16°37′S, 134°51′E) in the Northern Territory[34,35]. As individuals from these populations show no evidence of genomic admixture between subspecies[31,32] we hereafter refer to them as ‘parental’.

mtDNA comparison between subspecies 

Changes in the mitochondrial genome are required but not sufficient for inter-lineage mitonuclear incompatibility. Even a single substitution in the mtDNA can lead to reproductive isolation (e.g., [36]). To identify substitutions with the potential to contribute to incompatibilities, we sequenced whole mitochondrial genomes from each long-tailed finch subspecies and counted the number of fixed differences and associated amino acid changes within 13 protein coding genes between them. We assembled mitochondrial genomes using whole genome sequence (WGS) data from both long-tailed finch subspecies (P. a. acuticauda, N = 14; P. a. hecki N = 11) and their sister species the black-throated finch P. cincta atropygialis (N = 11; Table S1). We used bcftools (version 1.9)[37] consensus after calling and incorporating variants relative to the zebra finch reference genome (GCA_003957565.4) in each sample. Ambiguous positions were output with associated IUPAC codes[43]. We exclusively used WGS data derived from DNA extracted from muscle (rather than blood) to ensure that they were enriched for reads from mitochondrial DNA and not reads from the nuclear encoded mtDNA (NUMT) copy. The mitochondrial genome of the zebra finch reference was used as an additional outgroup to help assign substitutions to their lineage of origin. We used the MITOS2 WebServer[38] to extract specific sequences for each mtDNA encoded gene, including all 13 protein coding genes, 22 tRNA and 2 rRNA and the control region. Each gene was aligned using MAFFT (version 7) web-aligner[39] and we manually examined each alignment with Geneious prime (version 2023.1.2) to count and characterize the distribution of fixed differences between each of the three members of Poephila (Table S2). For each protein-coding gene, we calculated the average level of sequence divergence between taxa (D_XY) using the program DnaSP (version 6 [60]).

We next evaluated the extent of non-synonymous change between long-tailed finch subspecies for each protein-coding gene. We generated lineage-specific consensus sequences for each gene and used the vertebrate mitochondrial genetic code to translate these to amino acid sequences. We again used data from the black-throated finch and zebra finch as outgroups to assign amino acid changes to the long-tailed finch subspecies of origin. We classified each amino acid substitution by any associated change in biochemical properties. Finally, we calculated the ratio of the rate of non-synonymous to the rate of synonymous nucleotide substitutions (dN/dS) between long-tailed finch subspecies for each protein-coding gene using the model M0 implemented in the program Easy-CodeML (v1.41)[40].

Breeding

Embryos were produced by rotating pairs of long-tailed finches through 20 outdoor aviaries (4.1 m long x 1.85 m wide x 2.24 m high), with one pair per aviary containing nest boxes and nesting material. Mealworms, greens, dry seed, and water were provided ad libitum. We required females to have been exclusively with the partner for 14 days minimum before any experimental eggs were collected, to avoid embryos being sired by stored sperm from a different male. Nest boxes were checked every day and eggs were collected the day they were laid, given a unique number using a permanent marker, and stored in a ‘soft box’ at cool room temperature (22 ^oC) for up to seven days before being placed into an incubator (Brinsea Ovation 56 EX, Brinsea Products, Winscombe, U.K.). Eggs were incubated at 37.5^oC and 65% humidity for 12 days. On the 12^th day (i.e., the day before expected hatch) eggs were removed and measurements made. Working with embryos meant we eliminated the opportunity for the effects of extended parental care by parents of different genetic backgrounds.

Crossing Design

We produced backcrossed offspring by first breeding F₁ hybrid females with an unadmixed parental male of each subspecies, sequentially and in random order. Resulting offspring were defined as “mitoconcordant” backcrosses if their father’s lineage did match their mother’s mitotype and were defined as “mitodiscordant” backcrosses if father’s lineage did not match their mother’s mitotype (Figure 2). We note that in the literature these groups are often referred to as ‘maternal’ and ‘paternal’ backcrosses, respectively. Parental pairs of each subspecies were set up to breed at the same time and used as controls.

Mitodiscordant backcrossing exposes incompatibilities by producing offspring who inherit most of their nuclear genome from a different lineage to that of their mitochondrial genome. This is particularly pronounced in the case of the long-tailed finch, as >98.5% of fixed nuclear differences between subspecies are on the Z-chromosome[41]. In our cross design, mitodiscordant backcross offspring inherit Z chromosomes that have had no opportunity for recombination between subspecies. If male, they receive one Z chromosome from their hemizygous F₁ mother and a Z chromosome of matching subspecies identity from their father. Females have only a single Z chromosome inherited from their father. Mitodiscordant backcrosses have mismatched mitochondrial and nuclear genomes, but also have a low level of admixture on autosomes in their nuclear genome, which would otherwise be a potentially confounding source of hybrid dysfunction. For this reason, we included mitoconcordant backcross individuals as controls, as their level of autosomal admixture is the same as mitodiscordant backcrosses, but their mitochondria and Z chromosomes come from the same subspecies (Figure 2). It is worth noting that these several different types of backcrossed hybrids are also highly relevant to the wild hybrid zone, where we previously described that ~68% (243 of 357) of individuals had some level of admixture[32].

The goal of our cross design was to create embryos in which the mitochondrial haplotype (i.e., mitotype) of the two subspecies was placed in three distinct contexts: 1) a nuclear genetic background exactly matching the mitotype lineage (parental crosses), 2) a nuclear genetic background mostly matching the mitotype lineage (mitoconcordant backcrosses), and 3) a nuclear genetic background mostly not matching the mitotype lineage (mitodiscordant backcrosses). The six resulting groups of embryos were thus 1. parental acuticauda (acuticauda mother  acuticaudafather), 2. parental hecki (hecki mother  hecki father), 3. mitoconcordant backcross acuticauda (F₁ mother with acuticauda mtDNA  acuticauda father), 4. mitoconcordant backcross hecki (F₁ mother with hecki mitotype  hecki father), 5. mitodiscordant backcross acuticauda (F₁ mother with acuticauda mtDNA  hecki father), and 6. mitodiscordant backcross hecki (F₁ mother with hecki mtDNA  acuticauda father).  Mothers used in parental crosses only produced a single clutch while F₁ hybrid mothers used for backcrossing usually produced two clutches: one with an acuticaudafather and one with a hecki father (however, 4 F₁ females only produced a clutch with one partner). All fathers were only used for a single clutch.

Heart rate measurement

After 12 days of incubation each egg was removed from the incubator and immediately placed on a Buddy digital egg monitor (Vetronic Services, Abbotskerswell, Devon, U.K.), to measure embryonic heart rate from within the egg (a proxy for whole body metabolic rate)[42]. A timer was started when the egg was removed from the incubator, and heart rate and corresponding time since removal were recorded at least twice over two minutes. We used these measurements to calculate the rate of cooling and the predicted heart rate halfway through the period (i.e., at 60 seconds). Five individuals where heart rate was only measured once were excluded from this dataset.

High resolution respirometry strategy

High resolution respirometry was performed with permeabilised tissue samples from these embryos using an Oroboros O2k (Oroboros Instruments, Innsbruck, Austria). High resolution respirometry allowed us to systematically stimulate, uncouple and inhibit certain combinations of mitochondrial complexes (Table 3; Figure S1). The biochemistry and genetic architecture underlying the core proteins of oxidative phosphorylation (OXPHOS) is highly conserved across eukaryotes and has been well described [17,85] with four protein complexes (CxI, CxII, CXIII and CxIV) establishing the proton gradient across the mitochondrial membrane, and CxV driving ATP synthesis (see also Figure 1). We measured the maximum capacity of OXPHOS system when fuelled by electrons either through complex I (CxI), or both CxI and CxII in tandem (OXP-I, OXP-I,II), and the maximum capacity of the ETS when accepting electrons from CxII, and both CxI and II in tandem (ETS-I,II, ETS,II). Comparing our experimental groups for these measures reveals whether CxI is compromised. We also assessed the maximum capacity for complex IV in isolation to consume O₂, a measure not reliant on CxIII, which can reveal if CxIII is compromised. We also quantified the component of respiration that is not available to fuel ATP synthesis (LEAK-I,II).

These measures enable us to test whether mitonuclear interactions impact mitochondrial aerobic metabolism in long-tailed finch hybrids. More specifically, they allow us to interrogate the source of any mitonuclear dysfunction, be it from the two complexes with amino acid differences between subspecies (CxI and III) or from translational issues related to fixed differences in the mtDNA encoded tRNAs or rRNAs [43-45].

High resolution respirometry technical steps

After the embryonic heart rate was measured, embryos were used for high resolution respirometry. Embryos were euthanised via decapitation and then weighed (Mettler Toledo, PB303-S/FACT, Colombus, Ohio, USA). There was no significant difference between experimental groups in embryonic survival (X²₅ = 4.1, p = 0.53; Figure S3B) or mass at day 12 of incubation (X²₅ = 7.3, p = 0.20; Table 5; Figure S3C). The head and body were then laterally cut in two, and the right-hand side was placed in 400 µl of phosphate buffered saline (PBS), homogenized with three strokes of a teflon-glass Potter-Elvehjem homogenizer (Wheaton, 5 mL), and stored at -80 for later analyses including DNA extraction for sex determination (described below). The left-hand side was weighed, suspended in a dilution of 100 mg/mL of MiR05 medium (0.5 mM Egtazic Acid (EGTA), 3 mM MgCl2, 50 mM K-lactobionate, 20 mM taurine, 10 mM KH2PO4, 20 mM Hepes, 110 mM sucrose, free fatty acid bovine serum albumin (1 g L-1), pH 7.1;[46,47]); and then homogenized and mechanically permeabilised using a teflon-glass Potter-Elvehjem homogenizer (Wheaton, 5 mL). All work was done on ice. This homogenate was centrifuged for 1 min at 100 g to pellet big tissue, and then 500ul of the supernatant was added to a chamber of an Oroboros O₂k high resolution respirometer set at 37°C. The respirometer had two chambers and we always ran two samples in parallel. After two samples were loaded and the system had equilibrated (oxygen consumption was stable), 2 µl of digitonin (10 mg/mL) was added to ensure cells were fully permeabilised. We injected 5 µl pyruvate (5mM final), 5 µl malate (2mM final), 10 µl glutamate (10mM final) as electron donors for CxI, and then measured oxidative phosphorylation when only CxI was receiving electrons (OXPHOS-I) by adding 5 µl adenosine diphosphate (ADP, 1.25 mM final). We next measured coupled respiration capacity through both CxI and II (OXPHOS-I,II) by adding 25 µl succinate (10 mM final). We then determined O₂ consumption associated with proton leak when electrons were being received by both CxI and II (LEAK-I,II) by adding 5 µl oligomycin (10 uM final). We next estimated the maximum capacity of the electron transport system (ETS-I,II) by titrating the mitochondrial uncoupler carbonyl cyanide m-chlorophenyl hydrazine (CCCP; at 2 mM concentration) by first one injection of 3 µl, followed by 1 µl steps. We next estimated the ETS capacity when only receiving electrons through CxII (ETS-II) by adding 1 µl rotenone (0.5 uM final). We then measured non-mitochondrial O2 consumption by injecting 1 µl of antimycin A (2.5 uM final) to inhibit CxIII. Finally, we estimated the respiration capacity of CxIV directly by first adding 5 µl TMPD (0.5 mM final) and 5 µl ascorbate (2 mM final), and then calculating the background chemical respiration rate by adding 50 µl azide (100 mM final). Across all runs the chambers were opened consistently at the same stages for reoxygenation. After each run 1 ml of the solution was taken from each respirometer chamber and frozen at -80 ^oC, and at the conclusion of all measures we performed a Pierce BCA protein quantification assay on this stored solution (ThermoFisher, Scientific, Waltham, MA, USA). Protein quantity was included as a covariate in statistical models (below) to account for slight differences in the amount of biological material in each sample. In most cases a different embryo sample was run in each chamber, but in 40 cases the same embryo was run in duplicate to assess how repeatable these measurements were (Table 3). When two different samples were run, the first was sitting in the respirometer for some time before the protocol was begun, and to account for any variability introduced by this we included a covariate in our models that was the length of time between when the embryo sample entered the chamber and when the measure was made.

Molecular sexing

We extracted DNA from the homogenized right-side of all embryos using a Gentra PureGene kit (Qiagen, Valencia, CA, USA) following manufacturer instructions. We sexed embryos by amplifying an intronic portion of the CHD gene that differs in size between the Z-linked and W-linked copy following published protocols[32,48]. Gel electrophoresis of the PCR product revealed the sex of the embryo as either male (one band, homozygous, ZZ) or female (two bands, heterozygous, ZW).

mtDNA copy number

We estimated mitochondrial density based upon mtDNA copy number for all samples using real-time quantitative PCR (qPCR) analysis following protocols established for passerine birds (Stier et al. 2019, 2022). While the amount of mtDNA in a mitochondrion can vary, mtDNA copy number is a reliable proxy of mitochondrial density. We estimated the relative mtDNA copy number by measuring the amount of mtDNA relative to the nuDNA for each sample by qPCR on a CFX96 Real-Time System (BIO-RAD). We used RAG1 as a representative single-copy nuclear gene (F: GCAGATGAACTGGAGGCTATAA, R: CAGCTGAGAAACGTGTTGATTC) and used COX1 as our representative mitochondrial gene (F: TCCTAGCCAACTCCTCAC, R: CCTGCTAGGATTGCAAAT). A complete nuclear copy of the mitochondrial genome (i.e., a NUMT) exists in Estrildid finches; for example, in the zebra finch reference genome (GCA_003957565.4) at chr2:72232845-72249691. To minimize the risk of NUMT interference in estimating mtDNA copy number, we designed our COX1 primers to perfectly match the mitochondrial sequence observed in Poephila. Moreover, our use of embryo tissue is less likely to result in amplification of the NUMT as gDNA from muscle tissue has a much lower nuclear:mitochondria ratio than gDNA from blood[49]. We further verified that COX1 amplification resulted in a single product of the expected size on an agarose gel.

Following Stier et al. (2022), qPCR reactions had a total volume of 12 µl that included 6 ng of DNA, had PCR primers with concentrations of 200 nM and 6 µl of Absolute Blue qPCR Mix SYBR Green low ROX (Thermos Scientific). The qPCR conditions were: 15 mins at 95^oC, 40 cycles of 15 s at 95^oC, 30 s at 58^oC, 30 s at 72^oC. Samples were run in triplicate for each gene (RAG1 and COX1) across 5 plates, ten samples were pooled and included on every plate as a reference, and 20 samples were included on two separate plates. Amplicon efficiency was estimated from a standard curve of the reference between 1.5 to 24 ng (1.5, 3, 6, 12, 16, 24) run on one plate. The mean reaction efficiencies were 92.1% for RAG1 and 82.3% for COX1. To calculate relative mtDNA copy number for each sample we used the following formula: (1+Ef_COX1)^ΔCqCOX1/(1+Ef_RAG1)^ΔCqRAG1, with Ef as the amplicon efficiency, ΔCq as the difference in Cq between the reference and focal samples. Some samples were run twice on separate plates, with an inter-plate repeatability of 0.74 (95% C.I. [0.49, 0.88]), and based on triplicates intra-plate repeatability was 0.94 (95% C.I. [0.93, 0.95]). Seven samples were considered as outliers (Cq RAG1 < Cq COI) and removed from analysis.

Oxidative Damage assay

It is possible that some aspect of mitochondrial function other than aerobic respiration is impaired, but that the respiration rate is maintained. For this reason, we also measured the damage induced by reactive oxygen species (ROS). High ROS levels resulting from mitochondrial dysfunction can damage macromolecules, cell components and structures[43,44] and ROS-induced DNA damage correlates with shorter lifespan in the closely related zebra finch[45]. 8-OHdG (8-hydroxy-2’-deoxyguanosine) is a predominant form of ROS-induced oxidative lesions to genomic DNA that has been widely used as a biomarker for oxidative damage[50], and has been associated with increased mitochondrial respiration in birds [51,52]. While 8-OHdG does not block DNA replication, it contributes to mutagenesis and can affect DNA transcription[53]. Here we quantified levels of 8-OHdG in the DNA extracted from embryos using a competitive immunoassay (300 ng DNA, EpiQuick 8-OHdG DNA Damage Quantification Direct Kit Colorimetric, Epigentek, USA) following the manufacturer recommendations. This was run on a subset of samples (N = 182) chosen blind to respiration capacity to have roughly balanced numbers from each experimental group (N = 28-33 samples per group). The intra-plate coefficient of variation based on duplicates was 6.3 + 0.83. The inter-plate coefficient of variation based on four samples repeated over all three plates was 22.7 + 4.05, so we included plate ID as random effect in our analyses.

Statistical analyses

For statistical analyses we used R version 4.1.2 (R Core Team 2020) using RStudio version 1.2.5033 (RStudio Team, 2020) for the graphical interface. We used the lme4 and lmerTest packages[54,55] to run linear mixed models to test whether our measures of respiration capacity differed between the six experimental groups (parental acuticauda, parental hecki, mitoconcordant and mitodiscordant backcross with acuticauda mtDNA, mitoconcordant and mitodiscordant backcrosses with hecki mtDNA). These measures were each included as the dependent variable in separate linear mixed models, in which the fixed effects were experimental group, embryo sex, protein quantity in the chamber, and the length of time after embryo entered the oroboros chamber when the measure was made. In each we tested for an interaction between experimental group and embryo sex, and if this was not significant the interaction term was removed (but left in the model if significant). In each model parental pair ID was included as a random effect, to account for multiple members of the same clutch being included. Where there was a significant effect of experimental group, we ran Tukey post-hoc comparisons using the emmeans package[56] to assess specifically which groups differed from one another. Model assumptions were checked using the DHARMa package in R[57].

We similarly ran linear mixed models to test whether heart rate, oxidative damage and mitochondrial copy number differed between our experimental groups, but in these cases the fixed effects were experimental group, embryo sex, and for the heart rate model embryo mass, and the random effect was still parental pair ID. Oxidative damage and mitochondrial copy number were log-transformed to better meet model assumptions. An interaction between experimental group and sex was also included for each model, but then removed if not significant. In the mitochondrial copy number and oxidative damage models, plate ID were additionally included as random effects.  Embryo ID was also initially included as a random effect in the mitochondrial copy number and oxidative damage models, but it explained close to zero variance so was removed for the final models, likely because very few individuals were measured twice (10 for mitochondrial copy number, 21 for oxidative damage). A small number of values for mitochondrial copy number were extremely high (N = 10, total N = 211), well above the normally distributed other samples and were likely technical errors, and so were excluded to facilitate modelling.

References

17. Burton RS, Pereira RJ, Barreto FS. 2013 Cytonuclear Genomic Interactions and Hybrid Breakdown. Annu Rev Ecol Evol Syst 44, 281–302. (doi:10.1146/annurev-ecolsys-110512-135758)

31. Hooper DM, Griffith SC, Price TD. 2019 Sex chromosome inversions enforce reproductive isolation across an avian hybrid zone. Mol Ecol 28, 1246–1262. (doi:10.1111/mec.14874)

32. Lopez KA, McDiarmid CS, Griffith SC, Lovette IJ, Hooper DM. 2021 Evaluating evidence of mitonuclear incompatibilities with the sex chromosomes in an avian hybrid zone. Evolution 75, 1395–1414. (doi:10.1111/evo.14243)

36. Meiklejohn CD, Holmbeck MA, Siddiq MA, Abt DN, Rand DM, Montooth KL. 2013 An Incompatibility between a Mitochondrial tRNA and Its Nuclear-Encoded tRNA Synthetase Compromises Development and Fitness in Drosophila. Plos Genet 9, e1003238. (doi:10.1371/journal.pgen.1003238)

37. Li H. 2011 A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27, 2987–2993. (doi:10.1093/bioinformatics/btr509)

38. Donath A, Jühling F, Al-Arab M, Bernhart SH, Reinhardt F, Stadler PF, Middendorf M, Bernt M. 2019 Improved annotation of protein-coding genes boundaries in metazoan mitochondrial genomes. Nucleic Acids Res 47, 10543–10552. (doi:10.1093/nar/gkz833)

39. Katoh K, Rozewicki J, Yamada KD. 2017 MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief. Bioinform. 20, bbx108-. (doi:10.1093/bib/bbx108)

40. Gao F, Chen C, Arab DA, Du Z, He Y, Ho SYW. 2019 EasyCodeML: A visual tool for analysis of selection using CodeML. Ecol. Evol. 9, 3891–3898. (doi:10.1002/ece3.5015)

41. Singhal S, Leffer EM, Sannareddy K, Turner I, Venn O, Hooper DM, Strand AI. 2015 Stable recombination hotspots in birds. Science , 1–6.

42. Sheldon EL, McCowan LSC, McDiarmid CS, Griffith SC. 2018 Measuring the embryonic heart rate of wild birds: An opportunity to take the pulse on early development. The Auk 135, 71–82. (doi:10.1642/auk-17-111.1)

43. Martindale JL, Holbrook NJ. 2002 Cellular response to oxidative stress: Signaling for suicide and survival*. J. Cell. Physiol. 192, 1–15. (doi:10.1002/jcp.10119)

44. Birben E, Sahiner UM, Sackesen C, Erzurum S, Kalayci O. 2012 Oxidative Stress and Antioxidant Defense. World Allergy Organization J 5, 9–19. (doi:10.1097/wox.0b013e3182439613)

45. Marasco V, Stier A, Boner W, Griffiths K, Heidinger B, Monaghan P. 2017 Environmental conditions can modulate the links among oxidative stress, age, and longevity. Mech Ageing Dev 164, 100–107. (doi:10.1016/j.mad.2017.04.012)

46. Stier A, Romestaing C, Schull Q, Lefol E, Robin J, Roussel D, Bize P. 2017 How to measure mitochondrial function in birds using red blood cells: a case study in the king penguin and perspectives in ecology and evolution. Methods Ecol Evol 8, 1172–1182. (doi:10.1111/2041-210x.12724)

47. Ton R, Stier A, Cooper CE, Griffith SC. 2021 Effects of Heat Waves During Post-natal Development on Mitochondrial and Whole Body Physiology: An Experimental Study in Zebra Finches. Front Physiol 12, 661670. (doi:10.3389/fphys.2021.661670)

48. Lee JC-I et al. 2010 A novel strategy for avian species and gender identification using the CHD gene. Mol Cell Probe 24, 27–31. (doi:10.1016/j.mcp.2009.08.003)

49. Sorenson MD, Quinn TW. 1998 Numts: A Challenge for Avian Systematics and Population Biology. Auk115, 214–221. (doi:10.2307/4089130)

50. Valavanidis A, Vlachogianni T, FIOTAKIS C. 2009 8-hydroxy-2′ -deoxyguanosine (8-OHdG): A Critical Biomarker of Oxidative Stress and Carcinogenesis. J Environ Sci Heal Part C 27, 120–139. (doi:10.1080/10590500902885684)

51. Stier A, Schull Q, Bize P, Lefol E, Haussmann M, Roussel D, Robin J-P, Viblanc VA. 2019 Oxidative stress and mitochondrial responses to stress exposure suggest that king penguins are naturally equipped to resist stress. Sci Rep-uk 9, 8545. (doi:10.1038/s41598-019-44990-x)

52. Stier A, Monaghan P, Metcalfe NB. 2022 Experimental demonstration of prenatal programming of mitochondrial aerobic metabolism lasting until adulthood. Proc Royal Soc B 289, 20212679. (doi:10.1098/rspb.2021.2679)

53. Markkanen E. 2017 Not breathing is not an option: How to deal with oxidative DNA damage. DNA Repair59, 82–105. (doi:10.1016/j.dnarep.2017.09.007)

54. Bates D, Mächler M, Bolker B, Walker S. 2015 Fitting Linear Mixed-Effects Models Using lme4. J. Stat. Soft. 67, 1-47\8. (doi:10.18637/jss.v067.i01)

55. Kuznetsova A, Brockhoff P, Christensen R. 2017 lmerTest Package: Tests in Linear Mixed Effects Models. Journal of Statistical Software 82, 1–26. (doi:10.18637/jss.v082.i13)

56. Lenth R. 2020 emmeans: Estimated Maringal Means, aka Least-Squares Means. R package.

57. Hartig F. 2021 Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models [R package DHARMa version 0.4.4].

60. Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. 2017 DnaSP 6: DNA Sequence Polymorphism Analysis of Large Data Sets. Mol. Biol. Evol. 34, 3299–3302. (doi:10.1093/molbev/msx248)

85. Hill GE, Havird JC, Sloan DB, Burton RS, Greening C, Dowling DK. 2019 Assessing the fitness consequences of mitonuclear interactions in natural populations. Biol Rev 94, 1089–1104. (doi:10.1111/brv.12493)