Experimental evidence that thermal selection shapes mitochondrial genome evolution

Mitochondria are essential organelles, found within eukaryotic cells, which contain their own DNA. Mitochondrial DNA (mtDNA) has traditionally been used in population genetic and biogeographic studies as a maternally-inherited and evolutionary-neutral genetic marker. However, it is now clear that polymorphisms within the mtDNA sequence are routinely non-neutral, and furthermore several studies have suggested that such mtDNA polymorphisms are also sensitive to thermal selection. These observations led to the formulation of the “mitochondrial climatic adaptation” hypothesis, for which all published evidence to date is correlational. Here, we use laboratory-based experimental evolution in the fruit fly, Drosophila melanogaster, to test whether thermal selection can shift population frequencies of two mtDNA haplogroups whose natural frequencies exhibit clinal associations with latitude along the Australian east-coast. We present experimental evidence that the thermal regime in which the laboratory populations were maintained drove changes in haplogroup frequencies across generations. Our results strengthen the emerging view that intra-specific mtDNA variants are sensitive to selection, and suggest spatial distributions of mtDNA variants in natural populations of metazoans might reflect adaptation to climatic environments rather than within-population coalescence and diffusion of selectively-neutral haplotypes across populations.

suggests that polymorphisms that accumulate across mtDNA haplotypes found in different spatial locations have been shaped by selection to the prevailing climate. This hypothesis was termed the "mitochondrial climatic adaptation" hypothesis in 2017 24 .
Currently, this hypothesis remains contentious, primarily because the conclusions of previous studies have, by and large, been based on correlations between mutational patterns in the mtDNA sequence and climatic regions, many of which have proven difficult to replicate in other or larger datasets 25,26 . We therefore decided to apply experimental evolution to test the mitochondrial climatic adaptation hypothesis, by determining whether multigenerational exposure of replicated populations of fruit flies to different thermal conditions leads to consistent changes in the population frequencies of naturally-occurring mtDNA haplotypes.
In the wild, different locally-adapted populations routinely come into secondary contact and hybridize, which enables selection of novel mito-nuclear genotypes that might be better suited to a new or changing environment 27 . Such an evolutionary scenario is likely to have become increasingly common in the Anthropocene, wherein humans have rapidly altered both climatic conditions and levels of habitat connectivity 28 . We reproduced such a hybridization event under controlled laboratory conditions by interbreeding two subpopulations of D. melanogaster, each adapted to thermal environments at a different end of an established and well-studied latitudinal cline 29,30 . It is thought that the species was introduced into Australia during the past one to two hundred years, probably via recurrent introductions of flies from both African and European origins 30,31 . The species has been studied extensively in the context of thermal adaptation along latitudinal clines, both within Australia, and other replicated clines in other continents 29,30,32 . This research has shown that numerous phenotypes related to thermal tolerance exhibit linear associations with latitude, and that these patterns are underscored by linear associations of key candidate nuclear genes 29 . Yet, no research had focused on the quantitative spatial distribution of mtDNA variants 32 , until Camus et al. 24 reported that similar clinal patterns are found for two phylogenetic groups of mtDNA haplotypes (haplogroups) along the eastern coast of Australia. Furthermore, Camus et al. 24 were able to map these clinal patterns of mtDNA variation to the phenotype, showing that the mtDNA haplotype that predominates at subtropical latitudes confers superior resistance to extreme heat exposure, but inferior resistance to cold exposure than its temperate-predominant counterparts.
To characterize our model system in detail, we designed a study based on experimental evolution, in which we submitted replicated laboratory populations of D. melanogaster to one of four different regimes of thermal selection. We note that similarly to patterns observed in mtDNA haplotype frequencies, Wolbachia infection frequencies also concord to latitudinal clinal patterns along the Australian east coast distribution of D. melanogaster, with higher frequencies in low latitude populations 33 . Furthermore, both the mtDNA and Wolbachia are maternally-inherited. Therefore, it is possible that previously reported clinal patterns in mtDNA haplotypes in Australia 24 might have been in part shaped by direct selection on Wolbachia genomes, with changes in mtDNA haplotype frequencies brought about by genetic hitchhiking on particular strains of Wolbachia. In order to test the interacting effect of Wolbachia infection on the dynamics of mtDNA adaptation under thermal selection, we replicated our experiment, under two different conditions -one in which the ancestors of our experimental flies had been treated with antibiotics to remove Wolbachia infections, and the other in which the ancestors had not received antibiotic treatment.

Material and Methods
Experimental procedures. Wild subpopulations of D. melanogaster were sampled during March 2012 in Australia. We sampled a "hot" adapted subpopulation ("H"; Townsville: −19. 26, 146.79) in the northeast, and a "cool" adapted subpopulation ("C"; Melbourne: −37.99, 145.27) in the south of the continent. We collected fertilised females and established 20 isofemale lineages from each wild population. Each lineage then underwent three generations of acclimatisation to laboratory conditions. Wild fruit flies are often hosts of intracellular parasites, such as Wolbachia and associated maternally-transmitted microbiomes that are known to manipulate host phenotypes and affect their thermal sensitivity [34][35][36] . To assess the effects of thermal selection on the standing mitochondrial variation in our experiment, both in the presence and especially the absence of these maternally-inherited microbiota that co-transmit with the mtDNA, we treated a full copy of our isofemale lineages with the antibiotic tetracycline hydrochloride (0.164 mg mL −1 tetracycline in food for 3 generations), such that we maintained a one full copy with putative Wolbachia and unperturbed microbiomes, and one full copy without Wolbachia but with perturbed microbiomes 37 .
We then propagated these lineages for a further 10 generations to mitigate any effects of the antibiotic treatment under laboratory conditions. Flies were reared at 25 °C on a 12:12 hour light:dark cycle in 10 dram plastic vials on a potato-dextrose-agar medium, with live yeast added to each vial ad libitum. All isofemale lineages were then transferred from our laboratories in Australia to those in Japan, and their food medium changed to a corn flour-glucose-agar medium (see Supplementary Table S1), with live yeast added to each vial ad libitum. This transition to corn flour based food was made for logistic reasons and the macronutrient content of the two diets is qualitatively similar (Supplementary Table S1). We furthermore note that there was no scope for the isofemale lineages to have adapted to these potato-based food source provided in the laboratory in Australia, given these lineages were founded by single females, and thus lacked the requisite genetic variation required to facilitate an adaptive evolutionary response. Prior to setting up a series of replicated experimental populations, these lineages acclimatized to their new food source for a further 3 generations at 25 °C before entering the admixture process described below.
We  Table S2). We then allowed the flies to lay eggs over 8 consecutive days and transferred them to fresh bottles as indicated in Supplementary Table S3. From this step, we reared all flies in 250 ml bottles on a corn flour-glucose-agar medium until the experiment was concluded. In the next step, we created each experimental population by mixing F1 offspring (25 virgin males and 25 virgin females) from the HC bottles with corresponding F1 offspring (25 virgin males and 25 virgin females) from the CH bottles (25 ♂CH + 25 ♀HC + 25 ♀CH + 25 ♂HC). These flies (randomly selected from the first hybrid generation) we call the starting generation. In this way, we established 7 experimental populations from the tetracycline-treated isofemale lineages and 8 experimental populations from the untreated lineages. We allowed flies of these populations to mate and lay eggs at 25 °C before transferring bottles with approximately 500 eggs each into four thermal regimes, represented by cold versus warm temperatures on either a constant or fluctuating temperature cycle. This design nests the temperature treatment inside the experimental population; each experimental population consists of four experimental subpopulations, and each experimental subpopulation adapts to one of the four thermal conditions. Bottles maintained in the constant cold temperature were kept at 19 °C (denoted "19 °C"), and those in the constant warm temperature at 25 °C (denoted "25 °C"). In addition, we used Environmental Chambers (MIR-154, Sanyo) to generate fluctuating thermal conditions that are common in areas of origin of our experimental populations 38 Table S4). We controlled the size of each subpopulation by trimming egg numbers in each generation to approximately 500. At the end of the experimental evolution period of three months, adult flies were collected and fixed in 95% ethanol.
Data collection. Total genomic DNA was extracted using DNeasy Blood & Tissue Kit (Qiagen). To identify regions of variation between the populations, we sequenced total DNA of the H and C population samples quantitatively, using an Illumina platform at Micromon (Monash University, Australia). Length of reads was set to 200 bp, and we reached a minimum coverage of 26x at relevant sites of the mitochondrial genomes in each population 24 .
We mapped all reads to the published mito-genomic sequence NC 001709 in Geneious R6 39 . We observed overall mitogenomic variability and picked 14 mtDNA polymorphic sites (SNPs) that are not unique to the H or C populations. These SNPs segregate all flies into one of two corresponding mtDNA haplogroups denoted A and B in Camus et al. 24 (Fig. 2, Supplementary Table S5). Using multiplexPCR and MALDI/Tof mass spectrometry 40 (Geneworks, Australia), we genotyped individually a) all isofemale lineages; b) virtually all flies in the starting generation (50 males and 50 females in each of the 15 experimental populations) and c) at least 24 males and 25 females per experimental subpopulation in the final generation upon completion of experimental evolution. This method provides 100% accuracy in the assessment of individuals into one of the two haplogroups. A lower bound for the number of individuals per each sex and experimental subpopulation to sequence (n = 24) was estimated assuming α = 0.05 and a relative thermal effect of 10% at a power of 1 -β = 70%. In total, we genotyped 4410 individuals (data available in Supplementary Table S2). Data analysis. MtDNA is transmitted maternally; males do not transmit their mtDNA to their offspring; and therefore evolutionary changes in mitochondrial genomes must proceed via selection on females. Therefore, our analyses focus on estimating changes in mtDNA haplogroup frequencies in females of each experimental population across the three months of the experiment. We applied a linear mixed-effect model (lmer) in the lme4 v. 1.1-10 package 41 of R version 3.2.2 42 . in RStudio server v. 099.465 43 with restricted maximum likelihood estimation of variance components, and type III Wald F-tests with Kenward-Roger degrees of freedom appropriate for finite sample size 44 . We modelled the antibiotic and thermal treatment as fixed effects, and the experimental population (n = 15 levels) as a random effect. This model revealed a significant interaction between the fixed effects of antibiotic treatment and thermal regime on the frequency change of the B haplogroup (Table 1). To interpret the underlying basis of this interaction, we then decomposed the analysis into two linear mixed-effect models examining haplogroup frequency change according to thermal regime for descendants of flies treated by antibiotics (ATB) and untreated (UTR) separately (Table 2). We modelled the thermal treatment as a fixed effect, and the experimental population as a random effect. Statistical significance between fluctuating-cold and fluctuating-warm thermal treatments has been evaluated by multiple Welch's t-tests (Table 3). It is key to note that the four levels of the thermal treatment are not only characterised by differences in the environments (warm versus cold, in either constant or fluctuating states), but also by the opportunity for adaptation, given that the warmer treatments ran for seven generations, but the colder treatments ran for only three. While the effect of the number of generations of selection cannot be statistically partitioned from the effect of the thermal environments per se, one might plausibly predict that there will have been greater opportunity for changes in haplotype frequency under the warmer conditions, than under the cooler conditions. Importantly, however, we estimated selection coefficients for the thermal treatments, according to the haploid selection model 45 , that explicitly take into account generation time (Appendix S1). We also verified that our measurements are not affected by genetic drift by simulations of the Wright-Fisher model 46 (Appendix S1).
In order to evaluate whether haplogroup frequencies in males are following the frequencies in females in our starting generation, we applied a linear mixed-effect model comparison of haplogroup frequencies between males  and females in starting generation according to antibiotic treatment. Antibiotic treatment and sex were modelled as fixed effects; experimental population was modelled as a random effect.
Even within a single generation, thermal selection might yield sex differences in differential survival, from egg to adulthood, of individuals bearing different mtDNA haplotypes; however, under strict maternal inheritance these differences will be reset at each generation. To evaluate the capacity for thermal selection to evoke such within-generation sex-differences in mtDNA frequencies, we applied a multilevel model examining the effect of sex, antibiotic treatment, and thermal regime on B mtDNA haplogroup frequency in the final generation. We modelled the thermal treatment as a fixed effect, the experimental population (n = 15 levels) and subpopulations (n = 60 levels) as random effects. Statistical significance of within generation frequency differences between sexes in particular treatments has been evaluated by homoscedastic two-tail t-tests in Microsoft Excel.

Results
The A haplogroup is found to predominate in the low-latitude, hot, tropical subpopulation from Townsville (H), whilst the B haplogroup predominates in the temperate, cooler Melbourne subpopulation (C; Fig. 2 Table S2).
We observed a statistically significant two-way interaction between thermal regime and antibiotic treatment on changes in haplogroup frequency in our experiment ( Table 1). The interaction was driven by an effect of thermal regime on haplotype frequencies in the antibiotic treated, but not the untreated, populations (Group ATB, P = 0.0110, Fig. 3, Table 2). This effect is important because in the absence of Wolbachia infection, changes in haplogroup frequencies can presumably be attributed directly to selection on standing variation in the mitochondrial genome. In the antibiotic treated group, we found that the frequency of the B haplogroup decreased in both of the warmer treatments but increased in the colder treatments. This response is consistent with the spatial distribution of the haplogroups along the Australian cline, where the B haplogroup predominates at temperate higher latitudes, while the A haplogroup predominates in subtropical low latitudes 24 . The largest statistically significant haplotype frequency differences are observed between fluctuating cold and fluctuating warm conditions (P = 0.0034 in Table 3). We estimated the selection coefficient of the B haplogroup for fluctuating warm conditions s w = −0.082 ± 0.026. We estimated the selection coefficient of the B haplogroup for fluctuating cold conditions s c = 0.085 ± 0.050 (Fig. 4, Appendix S1). Simulations confirm that our observations are unlikely to be accounted for solely by drift (P = 0.0013 in Appendix S1).
We also observed an effect of the thermal regime on within-generation sex differences in the frequency of the mtDNA haplogroups (Tables 4 and 5, Figs 5 and 6). These differences were already apparent in the starting generation (P = 0.0387 in Table 4). The patterns observed in the antibiotic-untreated populations at 25 °C at the commencement of the experiment (panel I in Figs 5 and 6) tracked closely those observed at the conclusion of the experiment seven generations later for populations maintained under the same conditions (UTR 25 °C, panel II in Figs 5 and 6). This replication of the sex-specificity of mtDNA frequencies is striking, supporting the finding that sex-differences in frequencies are not occurring randomly (P = 0.0321 in Table 5). Under these particular conditions, the frequency of the B haplogroup was higher in adult males than in adult females, suggesting differences in egg-to-adult survivorship of the two haplogroups (UTR 25 °C, panels I and II in

Discussion
Our experimental evidence demonstrated context-dependent shifts in population frequencies of two naturally-occurring mtDNA haplogroups under thermal selection. This result is broadly consistent with the hypothesis that spatial distributions of mtDNA haplotypes in natural populations might be in part shaped by thermal selection. In experimental populations whose ancestors' coevolved microbiomes, including Wolbachia infection, were disrupted by antibiotic treatment (ATB), we observed significant effects of thermal selection on mtDNA haplogroup frequencies, in patterns consistent with their corresponding frequencies in the wild. The advantage of these conditions (ATB) is that we could attribute frequency changes of mtDNA haplogroups directly to the variation in the thermal selection regime. Furthermore, our estimated selection coefficients would predict that in the absence of demographic or genetic processes other than thermal selection in shaping the haplogroup distributions (Fig. 4), the mtDNA haplogroups would reach reciprocal mitochondrial fixation in less than 10 years on both sides of the eastern Australian thermal cline. stage, all flies had been maintained in standard laboratory conditions (25 °C) for 16 generations (14 generations as isofemale lineages, 2 during the admixture process). We then divided each of these 15 biological replicates into 4 subpopulations, subjecting each subpopulation to one of four thermal treatments (19 °C, 25 °C, fluctuating cold, and fluctuating warm), with each experimental subpopulation containing around 500 individuals. On the left side of the figure, yellow text denotes sample sizes associated with each stage of the admixture process for flies whose ancestors had been exposed to antibiotic treatment (ATB), while grey text on the right corresponds with untreated flies (UTR). Nevertheless, haplogroup frequencies have not reached fixation in wild populations along the Australian latitudinal cline 24 . Indeed, twenty years prior to the study of Camus et al. 24 , Boussy et al. 47 reported that both the A and B haplogroups coexisted within most sampled populations along the Australian eastern seaboard. Haplogroup A possesses the HinfI restriction site used in their study (see Supplementary Table S5). Clearly the dynamics of selection in the wild will differ from those in a highly controlled laboratory experiment, and spatial and temporal environmental variation is likely to lead to genotype-by-environment interactions that might maintain these mitochondrial haplotypes in the wild notwithstanding their sensitivity to thermal selection. In agreement with this argument, our experiments found that the predicted patterns of the A haplogroup outcompeting the B haplogroup under warmer conditions were sustained only in experimental populations that had been antibiotic-treated, and had therefore experienced a microbiome perturbation, and were free from Wolbachia. We did not observe the predicted patterns in populations in which flies harboured their coevolved microsymbionts, including possible infections with Wolbachia, which might have for instance lead to Wolbachia-induced cytonuclear incompatibilities 48,49 and diverse fitness related effects 34,50 . Wolbachia infection transmission dynamics are known to depend on thermal conditions 51 , and thus thermal selection in these untreated populations would presumably be mediated through complex interactions between host nuclear genome, mtDNA haplotype, Wolbachia genome, microbiont genomes, which might have obscured any direct effects on the mtDNA frequencies. Furthermore, in the wild, there is a latitudinal cline in Wolbachia presence 50   prevalence is likely to be itself shaped by climatic selection. The low-latitude Australian subtropical populations exhibit higher levels of Wolbachia infection than higher latitude temperate populations 33,50 . It is possible that genetic polymorphisms within the mitochondrial genome will interact with polymorphisms spanning distinct Wolbachia strains to affect fitness outcomes of their hosts. The nature of these complex interactions between Wolbachia, mitochondrial, and host nuclear genome currently remains completely unexplored. In the absence of further research that elucidates the relative contributions of selection on Wolbachia versus mtDNA haplotypes in shaping the patterns of mitochondrial haplotypic variation observed in nature, it remains difficult to derive predictions as to how host-Wolbachia dynamics will affect changes in mito-genomic compositions of natural fruit fly populations 33,52 . Wolbachia clades are also known to exhibit habitat-specific fitness dynamics 53 , and it is possible that different Wolbachia, or other microsymbiont, strains are linked to the two different mtDNA haplogroups studied here, given that each co-transmit with the mtDNA in perfect association along the maternal lineage, and that the mtDNA frequencies in the antibiotic-free treatments hitchhiked on frequency changes involving these microsymbiotic assemblages, as is expected by theory, and has been observed previously 54,55 . We observed that haplotype frequencies in males did not necessarily track frequencies in females across the experimental treatments. These sex differences across the experimental treatments were complex, and involved changes in sign under the different combinations of experimental conditions. A key point to note is that selection on mtDNA in males will not directly contribute to shaping patterns of mtDNA variation between generations, under the assumption that males virtually never transmit their mtDNA haplotypes to their offspring. As such, mitochondrial genomes are predicted to evolve under a sex-specific selective sieve 56 , in which mutations in the mtDNA sequence that confer harm to males can nonetheless accumulate in wild populations, so long as these same mutations are neutral or beneficial for females [57][58][59][60] . In the absence of inter-sexual positive pleiotropy, such male-expression specific mtDNA mutations could in theory shape patterns of haplotype frequencies within a generation, if they affect male-specific patterns of juvenile or adult survival, but would not be passed on to the next generation, and would thus not shape haplotype frequencies across generations. Male-biased mitochondrial genetic effects on key life history phenotypes have been observed previously 56,60-64 , and the patterns observed here strengthen the emerging view of ubiquity of sex-specific effects of mtDNA polymorphisms and suggest that such polymorphisms are sensitive to selection imposed by the thermal and microbial environment.
It is possible that paternal leakage (transmission of mtDNA haplotypes from fathers to offspring) could have affected the evolutionary dynamics of selection across the experimental treatments. While paternal leakage has been previously observed in Drosophila, most reports have been limited to cases involving interspecific crosses between individuals of divergent species 65,66 , or intraspecific crosses between individuals with high levels of genetic divergence 67 . One recent study, which sampled flies from natural European and Mediterranean populations of D. melanogaster, noted that as many as 14% of individuals were heteroplasmic for divergent haplotypes, thus indicating paternal leakage 68 . Notwithstanding, the average frequency of the minor haplotype was very low within individuals (less than 1%), and such a low level of heteroplasmy seems unlikely to be of evolutionary significance. Indeed, the Dowling lab had not observed even a single case of paternal leakage among the numerous mitochondrial strains created or maintained in the laboratory over the past decade, despite continual backcrossing to males of isogenic strains possessing a different mtDNA haplotype 24,60,69 . However, in 2017, Wolff et al. 70 observed heteroplasmic individuals in two of 168 replicate populations (1%) following a large experimental evolution study in which flies of two divergent mtDNA haplotypes coexisted across 10 generations. In summary, while paternal leakage appears to occur in this species at low frequencies, we assume at this stage that selection on males would have played at most only a minor role in shaping the intergenerational changes in the frequencies of each haplotype.  Mitochondrial genetic markers remain an important tool for population genetics, despite growing experimental evidence that mitochondrial genetic variation is affected by thermal 24 , and other kinds of selection 71 . The evolutionary trajectories of distinct mitochondrial haplotypes might furthermore be selected together with functionally-linked nuclear gene complexes 72,73 . This reinforces the point that phylogenetic, population-genetic, and biogeographic studies involving mtDNA should incorporate statistical tests to investigate the forces shaping sequence variation and evolution 74 , and examine variation at multiple genetic loci 75 . To date, researchers have focused mainly on the effects of nonsynonymous mutations in the evolutionary dynamics of mitochondrial genomes 76 , but growing evidence suggests that mitochondrial molecular function is also affected by single nucleotides in synonymous and non-protein coding positions on mtDNA 24 ; a contention supported by the current study since there are no non-synonymous SNPs separating the A and B haplogroups 24 .
Our study advances understanding of the dynamics of evolutionary adaptation by providing experimental evidence that thermal selection acts upon standing variation in the mtDNA sequence. However, further research is needed to resolve the dynamics of this thermal evolution; for instance, by determining whether thermal selection acts on the mtDNA sequence directly, or on epistatic combinations of mitochondrial-nuclear genotype or mitochondrial-microbial genomes 77 . Furthermore, it remains unclear how much of the pool of   non-neutral genetic variation that delineates distinct mitochondrial haplotypes has actually been shaped by adaptive relative to non-adaptive processes. Finally, because of the difficulty of implementing experimental evolution in vertebrates, almost all experimental work investigating the adaptive capacity of the mitochondrial genome has been conducted on a small number of model invertebrate species 19,59,71,[78][79][80][81] , with few exceptions [82][83][84] . Future studies should involve a combination of ecological and experimental evolutionary approaches, with high resolution transcriptomics and proteomics applied more generally across eukaryotes, and the development of tests enabling us to reliably discern the footprint of thermal selection in wild populations 85 .
Impact Summary. We applied experimental laboratory evolution to provide the first direct test of the "mitochondrial climatic hypothesis, " which predicts that the variation of mitochondrial genomes across natural distributions of metazoans can be shaped by thermal selection. Our design is the first of its kind when it comes to inferring the role of thermal selection in shaping mtDNA frequencies in nature. We harness two naturally occurring mtDNA haplotypes of Drosophila melanogaster that segregate along the east coast of Australia. One of these haplotypes predominates at subtropical northern latitudes and the other in the temperate and cooler south of the country. We then compete these haplotypes against each other in replicated experimental fly populations submitted to one of four different thermal regimes, in either the presence or absence of infection by Wolbachia, a coevolved endosymbiont that also exhibits maternal transmission. We confirm that when evolving in the laboratory under warmer conditions, a haplotype naturally predominating in subtropical conditions outcompetes a haplotype that predominates at cooler Australian latitudes in the wild. We see this effect on haplotype frequencies in females in populations where latent Wolbachia infections had been purged.
Our results also suggest that sex-specificity of mtDNA effects, and co-occurrence of other maternally-inherited microbiotic entities -of which Wolbachia is just one example -are likely to shape the trajectories of mitochondrial genome evolution in the wild.