Mitochondria as environments for the nuclear genome in Drosophila: mitonuclear G×G×E.

Mitochondria evolved from a union of microbial cells belonging to distinct lineages that were likely anaerobic. The evolution of eukaryotes required a massive reorganization of the 2 genomes and eventual adaptation to aerobic environments. The nutrients and oxygen that sustain eukaryotic metabolism today are processed in mitochondria through coordinated expression of 37 mitochondrial genes and over 1000 nuclear genes. This puts mitochondria at the nexus of gene-by-gene (G×G) and gene-by-environment (G×E) interactions that sustain life. Here we use a Drosophila model of mitonuclear genetic interactions to explore the notion that mitochondria are environments for the nuclear genome, and vice versa. We construct factorial combinations of mtDNA and nuclear chromosomes to test for epistatic interactions (G×G), and expose these mitonuclear genotypes to altered dietary environments to examine G×E interactions. We use development time and genome-wide RNAseq analyses to assess the relative contributions of mtDNA, nuclear chromosomes, and environmental effects on these traits (mitonuclear G×G×E). We show that the nuclear transcriptional response to alternative mitochondrial "environments" (G×G) has significant overlap with the transcriptional response of mitonuclear genotypes to altered dietary environments. These analyses point to specific transcription factors (e.g., giant) that mediated these interactions, and identified coexpressed modules of genes that may account for the overlap in differentially expressed genes. Roughly 20% of the transcriptome includes G×G genes that are concordant with G×E genes, suggesting that mitonuclear interactions are part of an organism's environment. © The American Genetic Association. 2021.

Genotype by environment interactions (G×E or GEI) emerge when different genotypes have different phenotypic responses to an environmental variable (Schmalhausen 1949, Via and Lande 1985, Scheiner 1993, Schlichting and Pigliucci 1998). Gene by gene interaction (G×G) or epistasis is when the effect of one gene on a phenotype is dependent on, or interacts with, the action of a different gene (Phillips 2008, Sackton and Hartl 2016). Biomedical science tends to prioritize “master genes” that cause certain traits (Lewis 1992). At the same time, it is evident that most traits of evolutionary and medical importance are continuous in nature and efforts to identify the genetic basis of this variation has focused on an additive model where many genes of small effect contribute to phenotypes (Boyle et al. 2017). G×E and G×G effects are often relegated to second-order status that complicate the identification of the genetic basis of phenotypes (Frazer et al. 2009). Studies that have sought to examine the influence of epistasis or G×E have found that these interaction terms can be as important as individual genes (Huang et al. 2012) or can significantly improve the accuracy of genomic prediction approaches (Arnau-Soler et al. 2019). While quantifying G×E and G×G effects has a robust statistical framework, the mechanistic basis of these terms remains unclear. By their context-dependent nature, epistasis and G×E make it difficult to predict phenotypes. This lack of predictability is a substantial problem for personalized or precision medicine, as well as for advances in agriculture and the challenges of climate change (Mackay et al. 2009). Moreover, G×E and G×G effects are commonly considered distinct phenomena, driven by external environmental versus internal molecular factors, respectively. As gene networks are increasingly understood as key components of the genetic architecture of complex traits (Furlong 2013, Boyle et al. 2017), it seems possible that G×E and G×G effects may be associated with the same genetic loci.

To study genome-wide epistasis or environment-wide G×E effects, one needs to determine the phenotypes of all genes in the genome with every other gene in the genome (or with every imaginable environment). This is not tractable: for an organism with 20 000 genes, 20 0002 = 4 × 108 tests are needed, a conservative estimate when all genes have multiple alleles in nature. To model this problem, we have focused on a subset of all interactions by studying combinations of mitochondrial and nuclear (mitonuclear) genotypes (Rand et al. 2018, Rand and Mossman 2020). The mtDNA encodes 13 protein coding genes that are subunits of the electron transport chain, plus 22 tRNAs and 2 rRNAs that translate these 13 mRNAs inside mitochondria. The nucleus contributes over 1000 gene products to the mitochondrion, enabling oxidative phosphorylation and ATP production, plus a variety of other metabolic and signaling processes (beta-oxidation of fatty acids, amino acid recovery, redox signaling, among others [Raimundo 2014, Shadel and Horvath 2015]). Mitonuclear communication involves anterograde (nucleus-to-mitochondria) signals and retrograde (mitochondrion-to-nucleus) signals (Liu and Butow 2006). As the central organelle responsible for a wealth of metabolic processes related to nutrient and oxygen flux, the combined mitonuclear genome is a nexus for epistatic (G×G) and G×E interactions (Guarente 2008).

Here we ask the question: are mitochondria environments for the nuclear genome? Using the Drosophila model, we pair alternative mtDNAs with alternative nuclear genotypes and expose these mitonuclear genotypes to alternative environments in factorial designs. Previous studies have identified mitonuclear G×E and G×G×E effects using various designs (Zhu et al. 2014, Aw et al. 2016, Mossman et al. 2016, Nagarajan-Radha et al. 2019). While these kinds of experiments have used pairings of mitochondrial and nuclear genomes that may be unlikely in nature (e.g., between allopatric populations or different species), they do model admixture events that might also be rare in nature (long distance dispersal, hybrid zones between differentiated forms, or admixture events such as human-Neandertal or -Denisovan). In addition, the experiments provide a test of the mitochondrial coadaptation hypothesis. In general terms, both of these kinds of scenarios predict that the degree of phenotypic effect of a mtDNA scales with the degree of nucleotide substitution between different mtDNAs. MtDNAs that commonly pair with nuclear genes are expected to provide a “native” genetic environment, while “foreign” mtDNAs that have never encountered a nuclear genome could generate some kind of incompatibility or novel phenotypic effect (Blier et al. 2001, Rand et al. 2004, Dowling et al. 2008, Sloan et al. 2017, Havird et al. 2019). We present new transcriptomic analyses of mitonuclear genotypes to address the question of whether G×G is independent of, or overlapping with, G×E effects. We compare these new analyses, that use a dietary switch as an environmental variable, with our recent results examining the effect of mitonuclear interactions on transcription under altered hypoxic conditions or nutrient signaling pathways associated with the target of rapamycin (TOR) complex. These analyses explore the hypothesis that mitochondria function as integrators of the pathways regulating both the internal genetic, and external physical, environments (Rand et al. 2018).

Materials and Methods

Drosophila Strains

Experimental mitonuclear genotypes were constructed by combining distinct mtDNAs from different female lines with distinct nuclear chromosomal complements from male lines. Details on the construction of these lines have been reported elsewhere (Mossman et al. 2016). The complete set of strains involved 6 mtDNAs paired with 12 nuclear genomes (nucDNA) from the Drosophila Genetic Reference Panel (DGRP; (Mackay et al. 2012, Huang et al. 2014)). The mtDNAs were derived from 3 wild strains of D. melanogaster (Oregon R [as “OreR”], Zimbabwe 53 [as Zim53], and AustriaW132 [as “Aut”]) and 3 wild lines of D. simulans (siI, siII [as sm21], and siIII [as “ma1”]). The “ma1” mtDNA is derived from D. mauritiana haplotype “mau12” that differs from D. simulans siIII by one base pair as a consequence of historical introgression between D. simulans and D. mauritiana (Ballard 2000, Ballard 2000, Montooth et al. 2010). The nuclear genomes were DGRP lines 304, 313, 315, 358, 375, 517, 707, 712, 714, 765, 786, 820 (see Figure 1A). Initial chromosomal replacements were done with balancer chromosome extractions onto each cytoplasmic background, followed by several generations of male backcrossing from the original DGRP line resulting in 6 × 12 = 72 mitonuclear genotypes (see Figure 1A). Wolbachia was removed from each female line prior to mitonuclear genotype construction. The notation for these genotypes is mtDNA;nucDNA, for example OreR;315 or sm21;820. These genotypes were scored for development time on 4 different diets (Mossman et al. 2016). We also report results from a different pair of mitonuclear genotypes where the OreR and sm21 mtDNAs are each placed on the OreR nuclear chromosomes using a similar balancer cross and backcross scheme OreR;OreR and sm21;OreR (Montooth et al. 2010, Villa-Cuesta et al. 2014, Santiago et al. 2021).

Figure 1.

(A) Experimental design for the construction of mitonuclear genotypes from different mtDNAs paired with specific DGRP nuclear genotypes. Three mtDNAs from geographic samples of D. melanogaster (OreR, Zim53, Aut) and 3 mtDNAs from the 3 distinct mtDNA haplotypes in D. simulans (siI, siII, siIII represented by “ma1”) were placed on each of 12 DGRP nuclear genomic backgrounds. The top 3 rows of the grid represent “Home Team” mitonuclear combinations, of D. melanogaster mtDNAs on D. melanogaster nuclear backgrounds (blue mtDNAs with blue nuclear chromosomes); the bottom 3 rows represent “Away Team” mitonuclear combinations of foreign D. simulans mtDNAs on D. melanogaster nuclear backgrounds (red mtDNAs with blue chromosomes). The phylogeny at the left shows the number of amino acid changes between the various lineages of the phylogeny. A prediction of the mitonuclear coadaptation hypothesis is that the Away Team combinations would show greater dysfunction due to mitonuclear incompatibilities resulting from ~2.5 million years of independent mutation in the respective genomes. (B) Experimental treatment showing the subset of mitonuclear genotypes used in the diet exposure and RNAseq experiment. The 4 mitonuclear genotypes developed from egg to adult on Lab Food. At age 4 days, replicate cultures of 30 single sex adults were transferred to fresh Lab Food for 2 days. At time = 0, replicates were frozen as controls, and remaining replicates were transferred to each of 2 alternative diets for 60 and for 120 min: Y = High P:C and to S = Low P:C for Yeast and Sugar respectively, and then frozen for later RNA extraction.

Diet Environments

The standard Lab Food contains 1.8 g agar, 5 g SAF yeast, 10.4 g yellow cornmeal, 22 g sucrose, 0.9 g tegosept in 4.5 ml 95% ethanol, cooked in 200 mL of distilled water (larger batches scaled up these ratios accordingly). Isocaloric diets were prepared from reciprocal mixtures of yeast extract and sucrose in a standard Drosophila food medium. The High Protein:Low Carbohydrate (High P:C), Equal P:C, and Low P:C diets have 32, 20, and 8 g of SAF yeast, and 3, 8, 20, and 32 g of sucrose, respectively (Mossman et al. 2016). All remaining ingredients are the same in each diet (1 g agar, 9 g yellow cornmeal, 0.45 g tegosept in 4.5 mL of 95% ethanol, cooked in 200 mL of distilled water). Development time was scored with development from egg to adult on each diet of these 4 diets (see Figure 2). The transcriptome analyses were conducted on a specific subset of mitonuclear genotypes raised on standard Lab Food to age 3- to 4-day-old adults, and switched to the High P:C or Low P:C diets for defined time periods of 0 hours (Control), 60 and 120 minutes as described below.

Figure 2.

G×G×E for development time in mitonuclear genotypes. The upper left depicts the 6×12, mito × nuclear genotypes placed on 4 different diets. The larger graph plots the development time for each mitonuclear genotype x diet combination, color coded by diet type. The insets show development times for mtDNAs on 2 specific nuclear backgrounds (DGRP-315 and DGRP-714) on each of the 4 diets.

Diet-Switch Experiments

Results are presented from 2 experiments: a new RNAseq analysis of a diet-switch experiment, and a specific RNAseq finding we reported in a recent publication focusing on mtDNA effects on the inhibition of TOR signaling using rapamycin. For the diet-switch experiment, flies were allowed to develop on Lab Food at controlled density and when adults reached age 3–4 days posteclosion, males and females of each genotype were sorted into 3 or 4 replicate single sex vials of 30 flies and allowed to recover for 2 days on Lab Food. Replicate vials of flies for the diet-switch treatment were then transferred to either High P:C (labeled Y for high Yeast) or Low P:C (labeled S for high Sugar) diets for 60 or 120 min, and then flash frozen in liquid N2. Control flies (labeled C) were from replicate vials that that were frozen at time 0 without any exposure to novel food. Each vial of 30 adult flies of one sex was the unit of replication for genotype, diet exposure, and RNA extraction. Figure 1B illustrates the procedures for this RNAseq diet-switch experiment that was performed on 4 mitonuclear genotypes. The transcriptional changes that are due to mito × nuclear G×G effects are inferred from factorial contrasts of mtDNAs and nuclear genotypes, pooling across all diet environments and time points. The transcriptional changes that are due to mitonuclear genotype × environment effects (G×E) are inferred from factorial contrasts of the 4 mitonuclear genotypes on standard Lab Food versus one alternative isocaloric diet (e.g., either High or Low P:C diets).

The results of this diet-switch experiment were compared with those from a distinct time-course transcriptome analysis using the OreR;OreR and sm21;OreR mitonuclear genotypes. The goal of the latter experiment was to assess the effects of rapamycin, an inhibitor of the TOR protein complex that regulates nutrient signaling, on gene expression in a diet and mtDNA-dependent manner (Santiago et al. 2021). Flies were raised on normal Lab Food, and adult males were separated into replicate vials of 30 flies at 3–4 days posteclosion to be used in a time-course experiment of joint rapamycin and refeeding. The intention of using rapamycin was to examine the effect of inhibiting the TOR complex that regulates nutrient signaling. All flies were starved on agar-water food overnight and one set of replicate vials were frozen in liquid N2. The remaining replicates were split between a treatment of rapamycin (200 µM in EtOH) or control (EtOH vehicle only) agar for 30 min. After 30 min, all flies were moved to normal Lab Food and frozen in liquid N2 as cohorts at 1, 2 , and 4 h poststarvation. This placed alternative mtDNA genotypes on or off rapamycin followed by refeeding. Additional details of this experiment are reported in a separate publication (Santiago et al. 2021).

RNAseq Analyses

Total RNA was extracted from 30 flies from each replicate vial using the Qiagen RNAeasy kit. RNA sequencing libraries were prepared by GeneWiz (South Plainfield, NJ) with a target size of 300 bp, and single-end reads were generated on an IIllumina HiSeq 2500 instrument, resulting in >20 million reads per sample. Fastq files were examined for quality using FastQC and libraries were filtered using the FASTX toolkit (v2.6) Reads were mapped to the dm3 reference genome of D. melanogaster obtained from the UCSC Browser using TopHat (v2.0.12) and Bowtie2 (v2.2.3). BAM files were converted to SAM files using samtools (v0.1.19), and read counts for annotated genome features were obtained using htseq-count in the HTSeq program (Anders et al. 2015). Details of parameters for these packages were reported previously (Mossman et al. 2019).

Statistical Analyses of RNAseq Results

Quantification and statistical analyses of genotype and diet effects on transcript levels were done using edgeR (Robinson et al. 2010) based on read-count data obtained as described above. In some cases, a replicate was dropped from the analysis if the RNAseq data did not pass filter or if that sample was a clear outlier in initial principle component and multidimensional scaling analyses. Four mitonuclear genotype were examined on different diets in the High P:C/Low P:C experiment: Zim53;DGRP-315, sm21;DGRP-315, Zim53;DGRP-820, sm21;DGRP-820 (Figure 1B). To determine the effects of mitonuclear genotype on nuclear gene expression, a full model of the fit was performed using the dispersion parameters for Common, Trended, and Tagwise conditions sequentially. To identify transcripts that were sensitive to mito × nuclear epistatic interactions, edgeR contrasts of the Zim53 versus sm21 mtDNA in the DGRP-315 background were contrasted against the Zim53-sm21 contrast in the DGRP-820 background separately in each sex, and diet treatment. For each diet and sex combination, we evaluated 3 different statistical models using the filtered read-count table for the libraries that pertain to the specific sex and diet treatments (Models 1, 2, and 3; see Figure 3A). Each model produces lists of genes for each term (e.g., mtDNA, nuclear, and mtDNA × nuclear), which can be ranked by P-value. We are interested in the overlap of the interaction term of Model 1 (a G×G list) with that of Model 2 (a G×E list). In the current analyses, we compare these models for the control food versus one specific alternative food (e.g., control vs. High P:C, or Control vs. Low P:C). Model 1 uses mtDNA and nuclear DNA each as factors with 2 levels and averages across the libraries for those genotypes exposed to different diets. Model 2 treats the combined mitonuclear genotypes as a “Genotype” factor with 2 × 2 = 4 levels while the “Environment” term is a factor that partitions the libraries between 2 levels: those from the control Lab Food treatment, and those from an alternative diet (e.g., High or Low P:C food). In each of Models 1 and 2, the same data set is used, but the data are partitioned differently. The same comparison can be done using Model 3, where the gene list from the mtDNA × nuclear (G×G) interaction term is compared with the gene list from the 3-way mtDNA × nuclear × diet (G×G×E) interaction term (see Figure 3A). These “overlap tests” are done for each sex and diet combination, for both the Models 1 versus 2 and the Model 3 contrast of 2-way versus 3-way interaction terms. Thus, with 2 sexes, 2 diet contrasts and 2 types of model contrasts, 8 possible “overlap tests” could be performed (see Supplementary Table S1). Each of these analyses identified nuclear transcripts with significant 2-way and 3-way interactions terms.

Figure 3.

Transcripts affected by mitonuclear G×G overlap with G×E transcripts.(A) Statistical models for assessing G×G and G×E effects on nuclear transcripts as a function of mtDNA (Mito), nuclearDNA (Nuclear), and Environment (Enviro). The interaction terms in red represent the focal statistic for comparing G×G and G×E effects. (B) Heatmap of coexpressed transcripts: rows are transcripts, columns are nuclear, mtDNA and diet and time treatments, as specified in the column headers. The C, Y, and S columns refer to Control Lab Food (C), or food with high Yeast or high Sugar, equivalent to High P:C or Low P:C food, respectively. The Time column headers refer to Control Food (0), or 1 (1) or 2 (2) h on the alternative food. (C) Venn diagrams showing overlap of G×G and G×E for females on the High P:C diet (top) or Low P:C diet (bottom), relative to lab food. (D) Permutation tests showing that observed data show greater enrichment of overlap (colored curves) compared to a randomization of the G×G and G×E gene lists. Data from 4 data partitions are shown: Females on High P:C and Low P:C and Males on High P:C and Low P:C diets, compared to Lab Food. The Venn diagram for males is in Supplementary Figure S1.

To test the null hypothesis that the G×G list is independent of the G×E list, we performed G-tests of heterogeneity using a 5% value for the null overlap, compared with the observed overlap of genes on these 2 lists. A 10% null was also used in addition, recognizing that equating a P-value of 5% in a typical statistical test may not equate with a 5% overlap of gene lists in the current context. In addition, we compared the G×G list to the G×G×E list from the full 3-way model to test the null overlap model. Details of this approach are described in (Rand et al. 2018) for a distinct experiment using hypoxia rather than diet as the environmental variable.

For the time course of the 2 OreR mitonuclear genotypes, we sought to identify transcripts that had distinct time-course profiles in the 2 genotypes exposed to the rapamycin versus control food. This was done using the R package impulseDE2 (Fischer et al. 2018) available in the Bioconductor suite of packages. Details of these analyses are reported in Santiago et al. (2021).

Functional annotation of the differentially expressed (DE) genes was inferred using the Gene Ontology enRIchment anaLysis and visuaLizAtion tool (GOrilla) package (Eden et al. 2009). Identification of coexpressed modules of genes was done using weighted gene coexpression network analysis (WGCNA) (Langfelder and Horvath 2008) using all libraries in the experiment detailed in Figure 3. To test for associations between the genes identified in coexpressed modules from WGCNA and those identified as “overlap” genes, we used a hypergeometric test to ask whether the identity of these genes is expected by chance as draws from all genes identified in the RNAseq analyses.

Results

Diet Environment Alters Mitonuclear Epistasis

Previously, we reported the development times for 72 mitonuclear genotypes built from 6 mtDNAs each placed on 12 different nuclear genomes of the DGRP. Three mtDNAs are from D. melanogaster (Dmel) and 3 mtDNAs from D. simulans (Dsim). The Dmel mtDNAs are OreR, Zimbabwe, and Aut, differing by ~8–10 amino acid changes across the ~3700 codons of the protein coding genes. The Dsim mtDNAs are siI, siII=sm21, siIII=ma1, differing by 10–40 amino acid positions. The divergence between Dmel and Dsim mtDNAs is ~100 amino acid positions (Ballard 2000). We have quantified development time for these 72 genotypes on 3 isocaloric Protein:Carbohydrate (P:C) diets, plus standard lab food (with lower P and C levels). This design allows direct partitioning of mtDNA, nuclear, and diet main effects, “mitonuclear” interaction (G×G) effects, genotype × environment (G×E) effects, and environmental modification of epistatic interaction (G×G×E) effects. Figure 2, showing all 72 genotypes × 4 diets = 288 combinations, reveals considerable G×E and G×G effects of this experiment. A high P:C diet accelerates development (Figure 2: yellow squares and lines), and low P:C diets delay development (Figure 2: red and white squares or red and black lines), but there are considerable changes of rank-order of the 72 genotypes across diets. These “crossing reaction norms” demonstrate G×G×E interactions, which can contribute as much to phenotypic variation as main effects (Mossman et al. 2016).

The line plot insets in Figure 2 show development time changes across 6 mtDNAs in 2 nuclear backgrounds (DGRP-315 and DGRP-714) and 4 diets. Comparable changes in the mtDNA “reaction norm” in each diet and nuclear background is evident for other nuclear backgrounds, and is supported statistically by a significant Mito × Nuclear × Diet (G×G×E) term in the ANOVA (Mossman et al. 2016). The additional insight provided by these figures—that yellow lines are “flatter” than black lines and that line patterns change between nuclear backgrounds—is that mitonuclear epistatic interactions are modified by diet.

Nuclear Transcripts Show Overlapping Mitonuclear G×G and G×E Effects

The pattern of development time G×G×E effects presented in Figure 2 reveals that a mitonuclear epistatic interactions can be modified by diet environment. The question arises whether the mitonuclear genetic interactions that alter development time are related to the genetic basis of diet sensitivity. Given the role that mitochondria play in both catabolic and anabolic nutrient processing, we chose to asses this association at the transcriptional level. To test this hypothesis, we have performed transcriptional profiling of mitonuclear genotypes under different environmental conditions. In a previous report, we showed that the transcripts most influenced by joint mitonuclear genetic effects (G×G) had significant overlap with transcripts influenced by joint mitonuclear genotype × hypoxia environment (Rand et al. 2018). We sought to verify this result by conducting an independent experiment with different mitonuclear genotypes, and a distinct environmental variable (diet, rather than hypoxia).

Figure 1B illustrates the design of the diet-switch RNAseq experiment. Figure 3B shows an example of a module of coexpressed genes that illustrate 3-way G×G×E effects and what common G×G and G×E genes can look like at the transcript level. In the DGRP-315 nuclear background, the 2 alternative mtDNAs (Dsim sm21 and Dmel Zim53) have no general differences, but each transcript is sensitive to diet (Figure 3B, heat map: rows are genes, Y and S at column headers specify high Yeast or high Sugar diets; High P:C and Low P:C, respectively). However, in the DGRP-820 nuclear background, the 2 alternative mtDNAs have dramatic effects on these transcripts: the sm21 mtDNA increases expression but is less responsive to diet, and Zim53 has reduced expression but is responsive to diet shifts.

If G×G and G×E effects are distinct processes, or are independent, these lists of transcript levels ranked by P-values should identify few overlapping genes among the highest-ranking transcripts (e.g., 5% by chance). If there are common factors to G×G and G×E, the overlap of these gene lists should be greater than random. The data reveal a significant enrichment for common genes: about 20% of transcripts are shared in these analyses: Figure 3C shows Venn diagrams based on the top 500 transcripts ranked by P-value for females experiencing a diet switch from Control food to either High P:C (Y) or Low P:C (S) food. Permutation tests reveal that the enrichment of these “overlap genes” is not an artifact compared to a randomized list (Figure 3D). To test the significance of this overlap, G-tests of heterogeneity were performed on the observed versus the 5% expected number of shared transcripts among the top-ranked genes, and observed overlap is significantly greater than that expected by chance (from Figure 3: Observed of 407:93 vs. a 5% expected overlap of 475:25: χ 2 = 43.13, P << 0.0001; using a 10% Expected of 450:50 gives χ 2 = 14.4, P < 0.00015; see Supplementary Table S1). Results of these overlap tests for males, and those using the Model 3 analyses are presented in Supplementary Figure S1. The results presented here are consistent with an earlier experiment with hypoxia. These 2 experiments used independent combinations of mtDNAs and nuclear genotypes and different environmental stressors (hypoxia or diet), and the ~20% “overlap” of G×G and G×E genes is repeatable (18–34% overlap for hypoxia; Rand et al. 2018).

Inferring Functional Associations of Overlap Genes

We examined the functional aspects of these “overlap genes” using the gene ontology (GO) platform GOrilla (Eden et al. 2009). GO analyses identified recurring themes that differed between the different statistical models and different food types. In the Model 1 versus 2 analyses (G×G vs. G×E), chitin and cuticle development was enriched as a “Process” and a “Function, while carbohydrate or lipid metabolism were enriched to varying degrees in either the High P:C or Low P:C diets. In the Model 3 analyses (G×G vs. G×G×E), immune or defense response and eggshell or vitelline membrane were common themes (see Figure 3 and Supplementary Figure S1 and Supplementary Tables S1 and S2). We note that the terms reported in Supplementary Table S1 have significant P-values, but some do not pass the false discovery rate threshold of 5%, and are provided for comparison only. Chitin metabolism, extracellular matrix, immune defense, and metabolism of lipids and carbohydrates are recurring terms in these analyses.

We further examined the association of these overlap genes with modules of coexpressed genes that were identified using WGCNA. The genes in the “steel blue” and “light cyan” modules identified in the WGCNA (see Supplementary Table S3) were significantly associated with the 93 overlap genes for females on High P:C (Y) food from models 1 vs. 2 (see Figure 3; P < 8.1e-07 for “lightcyan” and <4.1e-4 for “steelblue”). The genes in these 2 modules were also significantly associated with the 121 overlap genes in females on Low P:C (S) food (P < 1.2e-15 and < 2.6e-15, respectively, for these 2 modules; a Bonferroni threshold was set at 0.05/66 = 7.575e-4 based on tests of 66 WGCNA modules). These WGCNA modules are enriched for functions associated with carbohydrate metabolism and egg development, respectively (Supplementary Table S2); the “steel blue” module is shown in Figure 3B.

Mitochondrial Genotype Alters the Impact of TOR Signaling and Diet on Nuclear Gene Expression

The apparent overlap of G×G and G×E effects on nuclear expression under altered diets raises the question of how these effects are mediated. The TOR pathway regulates a broad set of interconnected pathways affecting nutrient signaling, autophagy, mitochondrial biogenesis among many other cellular activities (Saxton and Sabatini, 2017) . The drug rapamycin targets the TOR protein altering its kinase activity and multiple downstream effects. We have used rapamycin to explore whether alternative mitonuclear genotypes could modify the effects of rapamycin on gene expression during dietary supplementation. Previously we showed that mtDNA genotypes alter the impact of rapamycin’s effect on mitochondrial respiration and metabolite profiles (Villa-Cuesta et al. 2014). More recently, we have examined the effects of alternative mtDNAs and rapamycin on the transcriptional response to refeeding after starvation (Santiago et al. 2021). We compared the “home team” Dmel OreR strain carrying its own mtDNA (OreR;OreR) to an “away team” mitonuclear genotype carrying the Dsim sm21 mtDNA on the OreR nuclear background (sim21;OreR). Flies were starved overnight, then allowed to feed on control food, or food containing rapamycin, for 4 time points (0, 1, 2, 4 h). RNAseq analyses uncovered coexpressed modules of genes where the sm21 mtDNA inverts the temporal pattern of gene expression in response to rapamycin (Figure 4A: X-axis shows timepoints, genotypes, and control vs. rapamycin diet; significance inferred with Bioconductor package ImpulseDE2 [Fischer et al. 2018]). The module of coexpressed genes shown in Figure 4A is enriched for oxidative phosphorylation, carbohydrate and nucleotide metabolism, and several other metabolic processes (Santiago et al. 2021). While the use of 2 mtDNAs on one nuclear background precludes an epistatic G×G analysis, these analyses, using a distinct environmental exposure (diet + rapamycin), support the hypothesis that mitochondrial genotype is a modifier of the nutrient signaling effects of the TOR pathway and further demonstrate a strong G×E effect for mitonuclear genotypes.

Figure 4.

MtDNA genotype alters the transcriptional response to refeeding in a rapamycin-dependent manner. (A) Modules of genes that show inverted time course of expression for the “away team” genotype sm21;OreR when exposed to rapamycin (compare Rapa time courses for OreR;OreR vs. sm21:OreR). The transcripts in this module are associated with oxidative phosphorylation, carbohydrate metabolism, and nucleotide metabolism as shown in the heat maps below the time-course plots. (B) Cartoon showing that 2 distinct mitonuclear expression experiments both point to the transcription factor giant as a shared regulator of gene expression. (C) Network diagram illustrating the association between the giant transcription factor (left, red node) and the Dref transcription factor (right, blue node) and the respective transcripts that are associated with these transcription factors. Purple nodes are transcripts where analyses infer transcriptional modification from both giant and Dref.

What Are the Mechanisms of Shared G×G and G×E for Transcription?

The impact of mitonuclear genotypes, diet, and altered TOR signaling on nuclear transcript levels suggests that mitochondrial genotype alters the cellular environment in a manner that feeds back to transcriptional regulation. To explore this possibility, we use transcription factor binding site enrichment analyses to identified transcription factors that may be shared across the nuclear transcripts that are modified by mtDNA, and its interaction with nuclear and dietary factors. This was done for 2 independent experiments reported in (Mossman et al. 2019) and (Santiago et al. 2021).

The sets of DE transcripts in the experiment reported in Figure 4 (Santiago et al. 2021) were examined for shared transcription factor binding sites using the R-package RCisTarget (Aibar et al. 2017). For the cluster of coexpressed genes in Figure 4, 2 of the most highly enriched transcription factors were Dref (DNA replication-related element factor), which has been associated with mTOR activity, and gt (giant) a transcription factor involved in early developmental patterning in the Drosophila embryo. Together, Dref and gt were associated with 185 of the 258 genes in the enriched KEGG pathways shared in this coexpressed module of genes (Figure 4B and C).

Using an independent set of 8 mitonuclear genotypes, we performed RNAseq to identify nuclear transcripts that were DE in alternative mtDNA, or joint mitonuclear, genotypes (Mossman et al. 2019). These DE transcripts were mapped onto protein–protein interaction networks, and “neighborhood genes” nearby in the networks were also identified. Using a conservative set of transcripts that were shared between males and females, and 2 approaches for identifying gene neighborhoods in these networks, we used the oPOSSUM 3.0 package (Kwon et al. 2012) to identify enriched transcription factor binding sites. This analysis identified giant as significantly enriched among the DE transcripts (Figure 4B). It is notable that 2 different experiments, using different mitonuclear genotypes and experimental treatments, and different analysis pipelines and software packages both identified a common transcription factor as a possible mediator of the differential gene expression. This suggests that giant may be a key player in retrograde communication between mitochondria and nucleus.

Discussion

In this study, we have sought to use nuclear–mitochondrial interactions as a general model for exploring the relationship between genotype by environment interactions (G×E) and gene by gene interactions (G×G) or epistasis. Genetically, it is straightforward to pair different mtDNAs with different nuclear backgrounds, allowing direct tests of mitonuclear epistatic effects on phenotypes. Likewise, it is straightforward to expose these constructed mitonuclear genotypes to different environmental factors to quantify G×E. The intersection of these approaches makes the analogy that mitochondria are environments for the nuclear genome, and may provide a focused approach to explore relationships between G×E and G×G. The results presented here include both new unpublished findings, and specific aspects of previously published studies from our lab. Together the results seek to address the working hypothesis that mitochondria are integrators of G×G interactions influencing the internal cellular environment, and G×E interactions from external physical factors influencing organismal performance.

Quantifying G×G×E Is Not Enough

Mossman et al. (2016) quantified the joint effects of mtDNA, nuclear genome, and diet on development time in D. melanogaster. A full ANOVA model of the 72 genotypes (6 mtDNA × 12 nuclear) on 4 diets showed that diet had the greatest impact by far, with nuclear genotype having a moderate effect and mtDNA having a lesser effect. All interaction terms were highly significant, including the 3-way mtDNA × Nuclear × Diet effect, which contributed a small proportion of the overall variance (see Mossman et al. 2016, Table 1; Rand and Mossman 2018, Figure 2). The interaction plots shown in Figure 2 of this article illustrate the nature of the G×G×E interactions. The siI mtDNA slows development in the DGRP-315 nuclear background, and the accelerating effects of higher Protein:Carbohydrate diets cannot remove this mtDNA effect. But in the DGRP-714 nuclear background, the siI mtDNA developmental delay seen in the Lab Food is eliminated by a High P:C diet. Thus, the mitonuclear epistatic effect on development is altered by the diet environment. This can be quantified as a general pattern, implying that G×G and G×E have some functional connection, but the mechanisms underlying these effects cannot be inferred by further partitioning of the variance. Comparisons of the inferred phenotypic effect of amino acid substitutions in the various mtDNAs used identified mutations at position 313 in the ND2 protein and position 65 in the ND5 protein as impactful (Mossman et al. 2016, Supplementary Figure S1). However, these 2 positions are fixed between the D. melanogaster and D. simulans mtDNAs, so they cannot be assigned as uniquely causal for the siI-specific G×G×E effect. Because the mtDNAs and nuclear genomes used in the experiment are assayed as whole haplotypes, quantifying the statistical significance or effect size of G×G×E effects only provides a pattern. We need additional analyses to get closer to mechanism underlying these complex interactions.

Using Transcriptomes to Probe Mitonuclear G×G and G×E Effects

The use of transcriptional variation to quantify joint G×E and G×G effects is attractive for several reasons. First, it provides a large number of quantitative traits in one analysis: a typical RNAseq experiment in Drosophila can identify >12 000 traits with considerable accuracy given enough replication (Huang et al. 2020). Second, the level of a transcript from a specific gene is a quantitative trait, influenced by both external environmental factors that induce expression and internal molecular factors such as polymorphisms that influence initiation, spicing, and degradation of transcripts or aspects of chromatin availability (Gilad et al. 2008, Lopez-Maury et al. 2008). Third, annotated databases provide information about a gene’s function, facilitating mechanistic interpretation of quantitative variation. Finally, the coordinated expression of nuclear and mitochondrial genes interrogates aspects of central metabolism. These are likely to have general influence on most traits that contribute to organismal fitness in both a physiological and evolutionary sense.

We present new RNAseq data showing that transcripts exhibiting genotype-specific responses to changes in the protein:carbohydrate (P:C) ratio (dietary G×E transcripts) are shared with transcripts showing mtDNA-dependent responses to changes in nuclear genetic background (mitonuclear epistatic [G×G] transcripts; Figures 3 and 4). A previous study using an independent panel of mitonuclear genotypes and changes in oxygen environments (normal vs. 6% O2) also showed a statistically significant enrichment for transcripts with overlapping G×G and G×E effects (Rand et al. 2018). That study used 2 different mtDNAs (OreR from D. melanogaster and siI from D. simulans), plus different nuclear backgrounds (OreR and Austria from D. melanogaster). The new diet-based experiment confirms the general patterns of G×G and G×E overlap with entirely different genotypes and environments. Two observations are notable about the results from these 2 experiments. First, the proportions of transcript in the overlap set are similar between experiments, roughly 15–25%. Second, the heat maps shown in Figure 3 of this article and in the hypoxia experiment (Rand and Mossman 2018, Figure 6) show some general similarities. In both cases, the mtDNAs have clear impact on nuclear gene expression on one nuclear background, but not in the other (a G×G effect), while across the treatment conditions, the transcripts show genotype-specific sensitivity to the environmental variable (diet: this study or oxygen: previous study). This suggests some support for the general hypothesis that mitochondria are integrators of signals from both the internal genetic or molecular environment and changes in the external environment.

The assumption of these tests is that G×E and G×G transcripts would be independent, and the expected overlap between lists of genes from these interaction effects would be ~5%. G-tests of independence or chi-square tests reject this null, and do so if a 10% overlap is taken as the null expectation (see Supplementary Table S1). This indicates that the pattern of overlap is not likely a statistical artifact. The structure of the these separate 2-way G×E and G×G models (Models 1 and 2, see Figure 3) has the condition that the main effect of one model is part of the error term of the other model, making it unlikely that top-ranked genes in one model would also be top ranked in the other.

To examine functional aspects of the genes in the overlap set, we performed GO analyses using the GOrilla web application to compare the overlap set with the genes detected in our analyses. GO analyses identified recurring themes that differed between the statistical models and alternative food types. In the Model 1 versus 2 analyses (G×G vs. G×E), chitin and cuticle development were enriched as a “Function,” while carbohydrate or lipid metabolism were enriched to varying degrees in either the High P:C or Low P:C diets. In the Model 3 analyses (G×G vs. G×G×E), immune or defense response and eggshell or vitelline membrane were common themes (see Supplementary Table S2). Despite these differences, links between chitin metabolism, extracellular matrix, immune defense, and metabolism of lipids and carbohydrates emerge from these analyses.

Chitin is the primary structural component of the insect cuticle, formed from epithelial extracellular matrix as a composite of chitin and a variety of proteins (Moussian et al. 2015, Pesch et al. 2016). As the key component of the exoskeleton, the cuticle serves as the primary barrier between the organism and the environment, and as a dynamic structure, it is important in wound healing and homeostasis (Galko and Krasnow 2004). Chitin metabolism has been implicated in regulation of immune and inflammatory responses (Lee et al. 2011, Pesch et al. 2016), which may be the source of the GO annotation linked to immune and defense functions. Chitin is a polysaccharide composed of glucose monomers with one hydroxyl group on each sugar replaced with an acetyl amine group. The different annotation ranking for chitin-related functions between the High versus Low P:C diets may reflect different availability of amino acids and sugars for remodeling cuticle composition. The chitinases that regulate the breakdown of chitin exist as a multigene family in Drosophila, several of which function inside mitochondria (Von Ohlen et al. 2012). Because mitochondria serve as central organelles in processing dietary carbohydrates, including amino acid recycling and fatty acid oxidation, the different mitonuclear genotypes may well have different flux for these pathways. In turn, the different mitonuclear genotypes may modulate the activities of chitinases, either through intramitochondrial states or altered availability of substrates for chitin turnover.

Other recurring GO annotations were egg development, fatty acid biosynthesis, and carbohydrate metabolism. These terms were more common in the treatments with high sugar content in the diets (Low P:C). Excess sugar would certainly promote expression of carbohydrate genes, and excess sugar in the diet is known to stimulate fatty acid biosynthesis (Musselman et al. 2011). That different mitonuclear genotypes would modify these processes seems likely, and that these processes would be further altered by shifts in diet seems straightforward. While the link between these functions and chitin metabolism is speculative, these observations point to testable hypotheses about how combined mitonuclear genotypes may be modulators of the external and internal environment for biochemical pathways.

Are Coexpressed Modules of Genes the Basis for G×G and G×E Overlap?

Our statistical analyses of overlap between gene lists from G×G and G×E analyses assume independence of transcriptional control for each transcript. As modules of coexpressed genes are a common feature of gene expression profiles (Langfelder and Horvath 2008, Ayroles et al. 2009, van Dam et al. 2018), the assumption of independence may not hold. We used WGCNA (Langfelder and Horvath 2008) to identify modules of coexpressed genes in our data set (Figure 3; Supplementary Table S3) and then examined the correlation between the lists of overlap genes and those genes represented in coexpressed modules. The overlap genes identified in Figure 3C are significantly enriched for genes in the “steel blue” (Figure 3B) and “light cyan” modules identified in WGCNA, which are enriched for carbohydrate metabolism and egg developmental functions, respectively. This suggests that part of the shared expression aspects of mitonuclear G×G and G×E may lie in the nature of gene coexpression networks that regulate functions critical to mitochondrial pathways. G×E effects for transcription have been documented across most of the DGRP genotypes, and coexpressed modules of genes are largely preserved in the 2 thermal environments used in that study (Huang et al. 2020). In our mitonuclear transcription data, the magnitude and direction of change for each transcript is not the same in the G×G and the G×E analyses, but similar sets of genes respond as a WGCNA module. This could provide an explanation, if not a mechanism, for the statistical pattern of overlapping expression we have confirmed in these experiments.

What Are the Mitonuclear Communicators?

If mitonuclear G×G and G×E have some shared nuclear expression signal, what are the molecules that act as interorganellar “transcription factors?” The genetic designs of orthogonal pairing of mtDNAs with nuclear genomes implies that 2 criteria must be met: 1) the process must be mediated in trans (mtDNA genes are not cis with any nuclear gene) and 2) the signals must include retrograde directionality (from mitochondria to nucleus). This does not preclude secondary cis-regulated expression networks, or simultaneous anterograde signals from nucleus to mitochondria, but these mitochondrion-specific conditions must be part of the mechanism. There is considerable interest in identifying “mitokines,” or mitochondrial signals that regulate cellular function (Merkwirth et al. 2016, Conte et al. 2020). Likely candidates for signaling factors are small molecules such as NAD+/NADH ratios, ADP/ATP ratios or reactive oxygen species levels that may be linked to chromatin accessibility through histone deacetylases (Merkwirth et al. 2016). However, studies that have identified candidate mitokines are conducted on models with one mtDNA genetic background, so our current experiments employing mtDNA “alleles” may be able to contribute novel progress in this area.

As illustrated in Figure 4, 2 independent studies from our lab, using different mitonuclear genotypes and experimental conditions both identified the transcription factor giant as a shared factor across genes that showed altered transcription levels by mtDNA genotype. While it is not surprising that a transcription factor would be found as a shared communicator, since both data sets involved RNAseq, this convergence on giant is striking given the very different genotypes, and treatment conditions used. Giant is known as a transcription factor influencing expression of gap genes in early developmental events (Eldon and Pirrotta 1991). However, it is expressed through adulthood in males and females (Graveley et al. 2011), but at lower levels. Giant does have genetic interactions with TOR, which has been associated with mitochondrial function, but the role of giant in mitochondrial functions specifically has not been explored. It remains possible that the shared G×G and G×E effects that point to giant reflect the latter’s role in transcriptional regulation for many genes, increasing the chances or random overlap (Lawhorn et al. 2018). Nevertheless, these associations have now provided candidate factors to examine the functional link between genes as environments for other genes.

In summary, the results reported here provide support for the general pattern that mitonuclear epistatic interactions intersect with genotype-specific responses to environmental cues (Rand et al. 2018). The new data presented seek to determine the generality of this pattern and identify possible mechanisms underlying the pattern. The examination of GO categories enriched among the genes in the set of overlapping transcripts, and evidence that co-expressed modules of genes are part of this gene set, provides a possible functional explanation for the statistical pattern of overlap between G×G and G×E genes. If a module of coexpressed genes shows altered expression in response to any form of stress, the overlap of G×G and G×E gene lists may be a functional necessity. Further evidence that a common transcription factor is part of this expression network motivates future experiments to help clarify these mechanisms. Details aside, mitochondrial metabolism is likely to be a nexus for integrating epistatic and environmental factors of “stressors.” We characterize a G×E gene as a gene that is altered by the external environment, and a G×G gene as one that is altered by the action of another gene in the genome, or the internal environment. The boundary between the external and internal environments is blurred when one considers gene expression because all external signals must be processed through some internal pathways. Identifying these shared pathways is a key goal of integrating genotype with phenotype. Over 70 years ago, Fisher (1941) sought to characterize the components of genetic and environmental factors that contribute to effect sizes of gene substitutions, and placed all factors other than additive genetic effects, into a broad “environmental” term. Dominance, epistatic interactions, the environment, representing both biotic and abiotic effects, are all modifiers of the expression of a gene according to this model. While epistasis and G×E are commonly taught as different phenomena, there are solid quantitative and functional reasons why they may be shared factors. Mitochondria really are an environment for the nucleus, which has probably been true for more than a billion years before Fisher published his model.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.