Mitochondrial-Nuclear Interactions Mediate Sex-Specific Transcriptional Profiles in Drosophila.

The assembly and function of mitochondria require coordinated expression from two distinct genomes, the mitochondrial DNA (mtDNA) and nuclear DNA (nDNA). Mutations in either genome can be a source of phenotypic variation, yet their coexpression has been largely overlooked as a source of variation, particularly in the emerging paradigm of mitochondrial replacement therapy. Here we tested how the transcriptome responds to mtDNA and nDNA variation, along with mitonuclear interactions (mtDNA × nDNA) in Drosophila melanogaster We used two mtDNA haplotypes that differ in a substantial number of single nucleotide polymorphisms, with >100 amino acid differences. We placed each haplotype on each of two D. melanogaster nuclear backgrounds and tested for transcription differences in both sexes. We found that large numbers of transcripts were differentially expressed between nuclear backgrounds, and that mtDNA type altered the expression of nDNA genes, suggesting a retrograde, trans effect of mitochondrial genotype. Females were generally more sensitive to genetic perturbation than males, and males demonstrated an asymmetrical effect of mtDNA in each nuclear background; mtDNA effects were nuclear-background specific. mtDNA-sensitive genes were not enriched in male- or female-limited expression space in either sex. Using a variety of differential expression analyses, we show the responses to mitonuclear covariation to be substantially different between the sexes, yet the mtDNA genes were consistently differentially expressed across nuclear backgrounds and sexes. Our results provide evidence that the main mtDNA effects can be consistent across nuclear backgrounds, but the interactions between mtDNA and nDNA can lead to sex-specific global transcript responses.

THE ancient symbiosis that led to current day mitochondria and their eukaryotic host was a major milestone for intergenomic communication (Sagan 1967; Martin and Muller 1998). During the following ∼2 billion years, the resulting eukaryotic cell consolidated the genetic information from two ancestral genomes to a reduced set of genes in the mitochondrial organelle [mitochondrial DNA (mtDNA)] and many nuclear-encoded genes that are required for mitochondrial function. At ∼16.5 kb in size, the mtDNA encodes 13 oxidative phosphorylation (OXPHOS) proteins, 22 transfer RNAs (tRNAs), and 2 ribosomal RNAs (Taanman 1999). The remaining ∼1200 proteins associated with mitochondria are encoded by the nucleus (Wallace 1999). Mitochondria therefore require coordinated expression of genes from both mitochondrial and nuclear genomes for efficient function (Rand 2001; Smeitink ). Crucially, the effects of mtDNA variation, nuclear DNA (nDNA) variation, and their interactions on gene expression are poorly understood, yet they are expected to play a large role in an organism’s response to environmental or cellular stress.

Mitochondria perform many functions in the cell other than ATP production, including maintaining homeostasis, regulating redox signaling, and apoptosis (Friedman and Nunnari 2014). Recent studies have shown that when mitochondrial function is perturbed, retrograde signaling to the nucleus occurs via the unfolded protein response (UPR) (Houtkooper ; Quiros ), however the role of underlying mitonuclear genetic variation in retrograde signaling has not been studied in the context of gene expression. Aberrations in mitochondrial function are known to cause a wide spectrum of disorders in humans (DiMauro and Schon 2003; Lin and Beal 2006), and mitochondrial diseases with an mtDNA component can be difficult to characterize due to the complexity of dual genomic organization and large numbers of intergenomic interactions. The relatively high mutation rate, lack of protective histones, maternal inheritance, reduced effective population size, and no recombination all contribute to a higher frequency of deleterious mutations in mtDNA (Gemmell ). These “natural” genetic variants segregate together in complete linkage disequilibrium in the form of mtDNA haplotypes. Importantly, both de novo (spontaneous) mutations and haplotype variation have been associated with disease phenotypes across a number of species (Roubertoux ; Dimauro and Davidzon 2005).

While specific mtDNA mutations—associated with haplotype variation—increase the risks of some diseases, e.g., neurodegenerative diseases and aging (Stewart and Chinnery 2015), they do not operate in isolation from the substantial nuclear-encoded contribution to mitochondrial phenotypes (Wallace 1999; Wong 2012). Deleterious mitochondrial mutations often demonstrate incomplete penetrance (Giordano ), suggesting nuclear variants may dampen or amplify their effects. An unknown proportion of phenotypes demonstrating missing heritability are likely to involve the interactions (epistases) between mtDNA and nDNA (Mossman ). Indeed, mitonuclear interactions have been largely overlooked as sources of phenotypic variation (Pesole ), even though recent studies have shown mitonuclear interactions to be pervasive and context dependent in a suite of life history traits (Hoekstra ; Mossman ).

mtDNA variation is hypothesized to affect males more severely than females due to an imbalance in the strength of purifying selection between the sexes (Frank and Hurst 1996). Because mtDNAs are maternally inherited, males are at a genetical dead end for mtDNA evolution. Mutations that have negligible or zero effects in females, yet severe effects in males, can rise to high frequencies in populations by genetic drift (Frank and Hurst 1996). Under this hypothesis, it is predicted that mtDNA variation should promote larger phenotypic effects in males than females. Initial support for this hypothesis was suggested by mitochondrial disease cases that appeared more severe in males than females (Wallace 1992; Bernes ; Casademont ). More recent studies provide support for the Frank and Hurst hypothesis in Drosophila aging (Camus ) and fertility (Yee ; but see Friberg and Dowling 2008), yet there is a robust absence of support in development time and viability (egg-to-adult survival) (Mossman ). In humans, there is also equivocal evidence for the role of mtDNA haplotypes in fertility-related traits in males (Ruiz-Pesini ; Pereira ; Mossman ). There is some evidence in Drosophila that sex-specific mtDNA effects on gene expression exist (Innocenti ), and that variation in male-limited gene expression is dominated by mtDNA variation. However, it is not known whether differential gene expression associated with mtDNA variation is a robust phenotype across nuclear genetic backgrounds.

Here, we use a Drosophila mitonuclear introgression model (Montooth ; Meiklejohn ; Villa-Cuesta ; Holmbeck ) to test the hypotheses that mtDNA variation, nDNA variation, and mtDNA × nDNA interactions impact gene expression in males and females. A main motivation of this experiment was to examine whether different mitonuclear genotypes have distinct gene expression profiles. Here, we focused on a 2 × 2 genotype table (sensu Roubertoux ) to elucidate whether mtDNA variation per se, nuclear variation per se, and mitochondrial × nuclear variation influences gene expression in a sex-specific manner.

We have previously shown that across several Drosophila Genetic Reference Panel nuclear backgrounds, the siI and OreR mtDNA haplotypes diverge in phenotypic values (development time, egg-to-adult viability) (Mossman ). We have also shown there are nucleotide substitutions between siI and OreR mtDNAs that have putative deleterious effect on protein function (Mossman ), motivating an examination of gene expression variation that is associated with mtDNA and nuclear variation. We hypothesized that mtDNA variation preferentially modifies the expression of mitochondria-associated OXPHOS genes from both mitochondrial and nuclear genomes. Furthermore, we wanted to characterize the cellular processes and functional categories of genes that are modified by mtDNA, nuclear, and mitonuclear variation. We further hypothesized that mitochondrial protein translational machinery (Jacobs and Turnbull 2005) would be largely influenced due to its critical role in protein synthesis and the key participation of mitochondrial ribosomes. Finally, we investigated whether genes differentially expressed by mtDNA variation were enriched in male-limited genes (e.g., those involved in testes and sperm-related proteins) to test the generality of the Frank and Hurst hypothesis in alternative nuclear backgrounds.

Materials and Methods

Mitonuclear panel

Experiments were performed on four genotypes whose nuclear genomes were introgressed with mtDNAs from two sources: (i) OregonR strain of Drosophila melanogaster and (ii) siI strain of D. simulans. The nuclear backgrounds were Oregon R (OreR) and Austria W132 (AutW132), which are both D. melanogaster nuclear types. The four genotypes were (mito; nuclear): (i) OreR; OreR, (ii) OreR; AutW132, (iii) siI; OreR, and (iv) siI; AutW132. Specifically, the introgressions were performed by balancer chromosome replacement, followed by backcrossing to the original stock to homogenize the nuclear background and remove nuclear variation that may have been retained in the chromosome replacement process. The full double balancer replacement scheme is described in Montooth ). Alignments between mtDNA coding regions reveal substantial sequence divergence between OreR (NC_001709) and siI (AF00835) mtDNA haplotypes: a pairwise amino acid divergence (103 amino acids) and a pairwise synonymous divergence (418 SNPs) (Montooth ).

All strains used in this study were cleared of Wolbachia using tetracycline, and Wolbachia negative status was confirmed by PCR (Montooth ).

RNA sequencing sample preparation

Prior to RNA extraction, flies were reared on standard laboratory food at 25° on a 12 hr light:12 hr dark cycle, and density-controlled for one generation (five males and five females per vial). After density control, eclosed flies were allowed to mate for 3 days on standard food and were then separated by sex at a density of 50 males or 50 females per vial. After 2 days of recovery from CO2 anesthesia, batches of flies were transferred into a 1.5-ml microcentrifuge tube and immediately flash frozen in liquid nitrogen.

RNA was extracted from 30, 5-day-old whole healthy flies in each genotype by sex treatment. We used whole flies to test whether sex-limited gene expression was associated with mtDNA and nDNA genetic effects, as found by Innocenti . In addition, a number of other studies of genetic variation for gene expression in Drosophila have used whole flies, allowing for more direct comparison with our results (Gibson ; Ayroles ; Huang 2012; Mackay ; Huang ). In males, there were three replicates per genotype and in females there were three replicates in all genotypes apart from the siI; AutW132, which had two replicates. We followed a modified RNA sequencing (RNA-seq) sample preparation protocol from the Gilad Laboratory (Chicago University; http://giladlab.uchicago.edu/data/RNASeq_v2%202.doc). Messenger RNA (mRNA) was first extracted, followed by RNA fragmentation, complementary DNA (cDNA) first strand synthesis, second strand synthesis, end repair, poly adenylation, adapter ligation, and PCR enrichment. Throughout, RNA and DNA were quantified using the Qubit Kits (RNA Broad Range, dsDNA Broad Range, and dsDNA High Sensitivity) with a Qubit 1.0 Fluorometer. All Qubit reagents were obtained from Molecular Probes (Eugene, OR). Following PCR enrichment, we size selected PCR products with size range of 334–500 bp using a Caliper LabChip XT (DNA 750 chip) (Caliper Life Sciences, Hopkington, MA).

Gene expression assays

Gene expression was assayed in both sexes in all four genotypes using Illumina RNA-seq (HiSeq; Illumina, San Diego) with a 50-bp, single-end protocol. Samples were processed at Brown University’s Genomics Core Facility using an Illumina HiSeq2000 platform.

RNA-seq data preprocessing

RNA-seq read quality was first assessed using the FastQC v0.11.5 program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We then filtered the reads using a FASTQ quality filter (fastq_quality_filter) with –q 20 and –p 80 flags, as implemented in the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/commandline.html#fastq_quality_filter_usage). These values correspond to removing reads with 80% bases having a quality score of <20. Adapters were clipped from the reads using the fastx_clipper tool in the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/commandline.html#fastx_clipper_usage) and FastQC was repeated. We used TopHat v2.0.12 (https://ccb.jhu.edu/software/tophat/index.shtml) (Trapnell ) with Bowtie2 v2.2.3 (http://bowtie-bio.sourceforge.net/index.shtml) to map the reads to the dm3 reference genome using the dm3flybase.gtf annotation file obtained from the University of California Santa Cruz Genome Browser (https://genome.ucsc.edu/). BAM files generated by TopHat were converted to SAM files using SAMtools (http://samtools.sourceforge.net/) (Li ), and reads mapping to specific genome features (genes) were counted using htseq-count (http://www-huber.embl.de/users/anders/HTSeq/doc/count.html) (Anders ). The read-count data table obtained via HTSeq was used for downstream analyses. All preprocessing was conducted using computational resources and services at the Center for Computation and Visualization, Brown University.

RNA-seq data analysis

Differential expression (DE) analyses were conducted using two methods: DESeq (Anders and Huber 2010) following the vignette workflow (http://bioconductor.org/packages/release/bioc/vignettes/DESeq/inst/doc/DESeq.pdf), and edgeR (Robinson ). We also used a clustering algorithm to categorize genes based on their between-genotype expression profiles (see below). The measurement accuracy of RNA-seq was validated using a subset of genes and quantitative PCR (qPCR) as described in detail in Supplemental Material, File S1 (see also Table S1; Figure S13).

Gene filtering

In total we analyzed data generated from 23 RNA-seq libraries (11 female, 12 male). Throughout the data analyses, we used both complete data sets (including all genes), and filtered data (excluding the lowest 40th percentile of the data sensu the recommendations in the DESeq vignette). The purpose of the independent filtering was to remove genes from the total data set that have little or no chance of showing significant evidence of DE, while simultaneously increasing the detection power with a similar false discovery rate (FDR) (Bourgon ). We specify in each Results section which data set (unfiltered or filtered) and which analysis program was used.

Gene clustering

We used four genotypes in total for this study and our primary aim was to detect genes that show expression profiles consistent with mtDNA effects, nDNA effects, and mtDNA × nucDNA (epistatic) effects. We rationalized that different genetic effects would produce different norms of reaction shapes (profiles) across genotypes. We aimed to cluster genes based on their expression across different genotypes. For example, for all genotypes (independent variable), we asked which genes segregated with a similar pattern, and which of those patterns correspond with nuclear effects, mtDNA effects, and their mtDNA × nDNA interactions (see Figure S1 for theoretical gene expression outlines).

Model-based clustering of gene expression profiles was performed using MBClusterSeq (Si ), with k = 20 clusters in both sexes. Negative binomial (NB) models were used to perform the hierarchical clustering and hybrid tree builds as implemented in the MBClusterSeq package. Clusters are described using interaction plots of individual genes in each cluster and the mean genotype values per cluster were calculated. The log-fold change is relative to the normalized gene expression across all treatments (genotypes) with each row of the log-fold change matrix having a zero sum. The clustering algorithm allocates genes to a group based on their gene expression profile (in this case, the shape of the gene expression–genotype relationship).

Gene ontology enrichment

Gene ontology (GO) enrichment analyses were performed on clustered genes using the Bioprofiling.de program (http://bioprofiling.de/) (Antonov 2011) and the GO Consortium database (http://geneontology.org/) with the default submission and “Drosophila melanogaster” organism selected. The outputs of the “ProfCom_GO” analyses (Antonov ) were filtered for GO categories that evidenced a Bonferroni-corrected P-value of <0.05 and these GO categories were used in the heat map construction. In the sex-specific gene expression analysis, we used GOrilla GO enrichment (Eden ) to investigate the GO processes that were enriched in the male- and female-biased gene sets.

Analyses of OXPHOS gene subset

In addition to the analyses of the global gene set, genes encoding OXPHOS-related proteins (complexes I, II, III, IV, and ATP synthase) were selected for closer examination because these are putative targets of mitochondrial variation. Both mtDNA and nuclear genes targeted to the mitochondrion were downloaded from the MitoDrome database (Sardiello ; D’Elia ). Arc diagrams were used to display the relationship between various genes in the OXPHOS pathway in D. melanogaster (Tripoli ) and the effect of alternative mtDNA haplotype in each nuclear background and in each sex. DE was judged by DESeq in each nuclear background and sex separately. For clarity, only genes with a P-value of <0.1 are shown and a gene was only required to be included in one sex or nuclear background to be present in the arc diagram. A total of 78 nuclear genes and 13 mtDNA genes were included in the initial screen, of which 53 genes were included in the visualization (11 mtDNA genes and 42 nuclear genes; P < 0.1 in at least one genetic background).

Data availability

Drosophila strains used in this study are available upon request. FASTQ files from the 23 RNA-seq libraries are available at the Sequence Read Archive (SRA) under project accession SRP082430.

Results

Between-sex gene expression correlations

We first investigated whether there were signatures of sex-specific gene expression as a basis for understanding how mtDNAs may or may not modify gene expression in a sex-specific manner. Figure 1, A–D, shows biplots of the four mitonuclear genotypes using pseudocounts of the unfiltered data set. The pseudocount transformation is log2 (read count +1), and was calculated using the base mean in DESeq. There are clear gene subsets in the data distributions that show sex-biased expression, particularly for male-biased genes (Figure 1). The numbers of genes that are common across genotypic intersections are described for each sex (females, Figure 1E; males, Figure 1F). The female- and male-biased genes (red and blue data in Figure 1, A–D, respectively) were intersected separately and those elements that were common to all four genotypes were subjected to GO-enrichment analysis. GO-enrichment analyses (Eden ) show male-biased genes are enriched for sperm-related processes (Figure 1F and Table 1). Genes with female-biased expression (shown as red in Figure 1) are associated with, among other processes, egg-related GO categories. Table 1 describes the male-specific enrichment of GO categories.

Figure 1

Sex-biased gene expression across four mitonuclear genotypes. Gene expression profiles of individual genes are shown for each genotype analyzed in this study: (A) OreR;OreR, (B) siI;OreR, (C) OreR;AutW132, and (D) siI;AutW132. Biplots show female gene expression on the abscissa with corresponding male gene expression values on the ordinal scale. Data highlighted in red and blue show female- and male-biased genes, respectively. Sex bias was determined as a log2-fold >2 difference between females and males. Data in black show no sex bias in expression. (E and F) Venn diagrams describe the number of genes that are intersected between genotypes for sex-biased expression (red and blue genes in A–D). (E) Female and (F) male intersections are shown. Generally, there were more intersected genes that demonstrated male-biased expression than female-biased expression. Genes at the four-genotype intersection were subject to GO analysis. Male-specific GO processes are described in Table 1.

Table 1

Male-biased gene expression

GO term	Description	No. genes	Enrichment	Adjusted P-value
GO:0032504	Multicellular organism reproduction	75	3.58	1.82e−20
GO:0000003	Reproduction	75	3.43	1.75e−19
GO:0046692	Sperm competition	17	6.34	5.37e−08
GO:0044706	Multi-multicellular organism process	17	6.12	9.30e−08
GO:0048232	Male gamete generation	33	2.89	1.45e−05
GO:0007283	Spermatogenesis	32	2.83	3.60e−05
GO:0003341	Cilium movement	12	6.04	4.64e−05
GO:0048515	Spermatid differentiation	10	5.6	1.47e−03

Significantly enriched GO categories are shown for the genes that are intersected between all four mitonuclear genotypes. P-values were adjusted using the Benjamini and Hochberg (1995) method.

Genotype signatures of gene expression

There are large differences between the sexes in genotype signatures of gene expression. Figure 2 shows multidimensional scaling (MDS) plots (Ritchie ) of each genotype in female and male data sets. The data are plotted in two dimensions and the distances between samples (libraries) approximately reflect the leading log2-fold change between the samples for the genes that distinguish those samples (Ritchie ). Data were filtered as the top 10,000, top 1000, top 100, and top 10 genes. The top genes are those with the largest SD in expression between samples. In females, the progressive filtering increases the genotype-genotype distance, whereas in males, the progressive filtering increased the distance between mitochondrial genotypes only in the Aut132 nuclear background (Figure 2).

Figure 2

Genotype signatures of transcript variation. MDS plots show progressive filtering of genes based on their expression differences. The distance between points in a plot approximately reflects the relatedness between libraries based on their transcript measures. (A–D) Female and (E–H) male profiles are shown. The top 10,000 (A, E), top 1000 (B, F), top 100 (C, G), and top 10 (D, H) most deviant genes are shown. Separate genotypes are color coded: siI;OreR, yellow; OreR;OreR, black; siI;AutW132, blue; OreR;AutW132, red. Mitonuclear genotypes were clearly distinguishable across all filtering levels in females. Only the most differentially expressed genes were able to distinguish mtDNA haplotypes in males.

mtDNA substitution effects

We next tested whether mtDNA substitution conferred DE of genes in each nuclear background and in both sexes, separately. We did this using DESeq (Anders and Huber 2010) on individual nuclear backgrounds for each sex (two tests per sex). The read-count data were modeled using NB distributions. The fold change and associated P-value were calculated for each gene and P-values were adjusted to account for multiple testing using the Benjamini–Hochberg method (Benjamini and Hochberg 1995), providing FDRs. In all cases the Oregon R mtDNA was the reference mtDNA background and effects of the siI haplotype are reported in the analyses. For nDNA results, the OregonR nDNA was the reference and the reported fold changes are for AutW132/OregonR.

In both females and males, mtDNA substitution conferred DE of genes from both mtDNA and nuclear genomes (Figure 3). In general, the response to mtDNA substitution was greater in females than males. Both sexes and nuclear backgrounds responded with distinct patterns. Among the most consistently modified genes were those encoded by mtDNA, with siI haplotypes mainly conferring downregulation in gene expression relative to the OregonR baseline. The majority of mtDNA OXPHOS protein-coding genes were downregulated (comparison = siI/OreR) in both males and females and in both nuclear backgrounds, with exceptions for mtDNA ND2, which was consistently upregulated across sexes and nuclear backgrounds.

Figure 3

Effects of mtDNA substitution on gene expression across nuclear backgrounds and sexes. Volcano plots describe the log2-fold change in expression of genes and their corresponding −log10 P-value, as determined by DESeq. Female genotypes are shown on the top panel: (A) Oregon R nDNA, and (B) AutW132 nDNA. Males are shown on the bottom panel: (C) Oregon R nDNA, and (D) AutW132 nDNA. Data in red are the mtDNA genes. There are nuclear background effects on mtDNA substitution and females generally showed more effects of mtDNA haplotype on nuclear gene expression. Horizontal dashed lines show P-value cut offs at equivalent P = 0.05. Vertical dashed lines show ±2 × fold up- or downregulation of a gene due to alternative mtDNAs with values shown for siI mtDNA, relative to the OregonR mtDNA. Outliers are not shown and females have a different magnitude of variation on the ordinal scale than males. Note the consistent red datum in the top right section of each plot (ND2 gene).

Volcano plots in Figure 3 describe the effects of mtDNA substitution. Shown are the results of DESeq analyses for the unfiltered data set. The results for the filtered data set (top 60% quantile) are shown in Figure S2. In the unfiltered- and filtered-gene analyses, estimation of size factors to normalize the counts to a common scale was performed on the unfiltered and filtered data, respectively.

We then intersected the genes that were significantly differentially expressed by mtDNA substitution in each nuclear background (×2) and in both sexes (×2). In this analysis, we did not require a specific direction of effect, e.g., up- or downregulated, just that the genes were present at a P-value threshold. A summary of the number of genes that are present, at various P-value thresholds, are shown in Figure S3. Those genes that are consistently differentially expressed (at P < 0.05) by mtDNA substitution across nuclear backgrounds and sexes are described in Table 2. Specifically, we report a conservative analysis in which all four nuclear background × sex treatments are intersected (A∩B∩C∩D)—where ∩ signifies intersection—and a more parsimonious analysis, in which a gene was only required to be included in a three-way intersection, e.g., including A∩B∩C, A∩B∩D, A∩C∩D, and B∩C∩D (Table 2). In the strict four-genotype intersection, only five genes were significantly differentially expressed by mtDNA variation in (A) female OregonR background, (B) female AutW132 background, (C) male OregonR background, and (D) male AutW132 backgrounds. These were three mtDNA genes (ND2, ATP6ase, and ND4L), an unannotated computed gene (), and (Ser4; a serine protease). High within-group variance that results from pooling across other experimental factors (e.g., nuclear backgrounds) is likely to reduce the sensitivity to detect first-order effects of mtDNA, so the five core genes influenced by mtDNA variation is likely a conservative estimate of the true number.

Table 2

Genes intersected in both males and females and in both nuclear backgrounds in response to mtDNA substitution

Gene name	Flybase ID	Location	Symbol	CG number
A∩B∩C∩D intersection
CG11966	FBgn0037645	3R: 8,950,246..8,963,037 (−)	CG11966	CG11966
Mitochondrial NADH-ubiquinone oxidoreductase chain 2	FBgn0013680	mitochondrion_genome:240..1,265 (+)	mt:ND2	CG34063
Mitochondrial ATPase subunit 6	FBgn0013672	mitochondrion_genome:4,062..4,736 (+)	mt:ATPase6	CG34073
Mitochondrial NADH-ubiquinone oxidoreductase chain 4L	FBgn0013683	mitochondrion_genome:9,545..9,835 (−)	mt:ND4L	CG34086
Jonah 25Bi	FBgn0020906	2L: 4,954,279..4,955,144 (−)	Jon25Bi	CG8867
A∩B∩C intersection
Cyp28d1	FBgn0031689	2L: 5,210,460..5,212,445 (+)	Cyp28d1	CG10833
—	FBgn0035300	3L: 1,971,638..1,974,285 (+)	CG1139	CG1139
—	FBgn0038068	3R: 12,636,540..12,637,976 (+)	CG11600	CG11600
—	FBgn0033774	2R: 12,762,113..12,763,825 (+)	CG12374	CG12374
Chemosensory protein B 38c	FBgn0032888	2L: 20,820,089..20,820,867 (−)	CheB38c	CG14405
Integrin βν subunit	FBgn0010395	2L: 21,053,033..21,058,044 (+)	Itgbn	CG1762
εTrypsin	FBgn0010425	2R: 11,345,237..11,346,066 (−)	εTry	CG18681
Bent	FBgn0005666	4: 724,400..776,474 (+)	Bt	CG32019
—	FBgn0032494	2L: 13,239,989..13,241,773 (+)	CG5945	CG5945
—	FBgn0046999	2R: 16,873,015..16,873,755 (+)	CG6429	CG6429
Esterase 6	FBgn0000592	3L: 12,188,818..12,190,705 (+)	Est-6	CG6917
Hemolectin	FBgn0029167	3L: 13,846,054..13,860,001 (+)	Hml	CG7002
Henna	FBgn0001208	3L: 7,760,453..7,763,166 (+)	Hn	CG7399
Immune-regulated catalase	FBgn0038465	3R: 17,002,875..17,007,208 (+)	Irc	CG8913
A∩B∩D intersection
Phosphoenolpyruvate carboxykinase	FBgn0003067	2R: 18,536,767..18,539,416 (+)	Pepck	CG17725
Ugt36Bc	FBgn0040260	2L: 16,799,025..16,801,584 (+)	Ugt36Bc	CG17932
Cyp6d5	FBgn0038194	3R: 14,029,228..14,031,483 (+)	Cyp6d5	CG3050
	FBgn0052023	3L: 8,803,602..8,804,182 (−)	CG32023	CG32023
Ankyrin 2	FBgn0261788	3L: 7,655,389..7,718,395 (−)	Ank2	CG42734a
	FBgn0053346	3R: 28,669,233..28,670,452 (+)	CG33346	CG33346
lectin-37Db	FBgn0053533	2L: 19,419,075..19,419,767 (+)	lectin-37Db	CG33533
	FBgn0053965	3L: 1,232,882..1,234,015 (+)	CG33965	CG33965
	FBgn0085256	2R: 11,248,126..11,248,610 (+)	CG34227	CG34227
	FBgn0038820	3R: 20,676,224..20,677,351 (−)	CG4000	CG4000
	FBgn0039474	3R: 27,020,238..27,021,378 (−)	CG6283	CG6283
	FBgn0037936	3R: 11,688,381..11,690,051 (−)	CG6908	CG6908
	FBgn0039670	3R: 29,587,237..29,588,104 (−)	CG7567	CG7567
	FBgn0039687	3R: 29,724,159..29,724,950 (−)	CG7593	CG7593
B∩C∩D intersection
PGRP-SC2	FBgn0043575	2R: 8,716,950..8,717,695 (+)	PGRP-SC2	CG14745
PGRP-SC1a	FBgn0043576	2R: 8,709,733..8,710,320 (+)	PGRP-SC1a	CG14746
Mitochondrial NADH-ubiquinone oxidoreductase chain 3	FBgn0013681	mitochondrion_genome: 5,608..5,961 (+)	mt:ND3	CG34076
Mitochondrial cytochrome b	FBgn0013678	mitochondrion_genome: 10,499..11,635 (+)	mt:Cyt-b	CG34090
Mitochondrial NADH-ubiquinone oxidoreductase chain 1	FBgn0013679	mitochondrion_genome: 11,721..12,659 (−)	mt:ND1	CG34092

The genes required an unadjusted P-value ≤0.05 to be included in the analysis. The four-genotype intersection along with three three-way intersections are shown: (A) female OregonR background, (B) female AutW132 background, (C) male OregonR background, and (D) male AutW132 backgrounds. There were no genes present in the A∩C∩D intersection. Genes in boldface font are mtDNA protein coding genes.

CG42734 overlaps with two additional computed genes (CGs): CG44195 and CG32373.

Mitochondrial OXPHOS genes and mtDNA variation

mtDNA variation conferred different effects on global gene expression between the sexes and nuclear backgrounds (Figure 3). We next focused on the nuclear and mtDNA genes of the OXPHOS pathway, since these are jointly encoded by mtDNA and nuclear genes and hypothesized to be more sensitive to mtDNA and nuclear covariation. We filtered our global DE data sets for the OXPHOS genes and graphed those genes that were differentially expressed (Figure 4). We divided the genes into their respective OXPHOS complexes.

Figure 4

OXPHOS-related genes are differentially expressed by mtDNA substitution. Arcdiagrams (Sanchez 2014) show the identities of OXPHOS genes whose expression levels are affected by mtDNA variation. Line thicknesses correspond with the significance level. Gene identifiers in red are those encoded by mtDNA. Male effects (black lines) and female effects (gray lines) are shown. Different OXPHOS complexes are differentiated by circle color: complex I, black; complex II, blue; complex III, purple; complex IV, green; ATP synthase (complex V), gray. The effects in the OregonR nuclear background are shown in the left plot; AutW132, on the right. The data shown are from the filtered data set. nucDNA, nDNA.

In males, there was a high degree of symmetry in the effects of mtDNA variation across the nuclear backgrounds (Figure 4). The differentially expressed genes (with P < 0.05) were exclusively targeted to mtDNA-encoded proteins within the various complexes (red gene identifiers in Figure 4) across both OregonR and AutW132 nuclear backgrounds. In the AutW132 nuclear background, fewer mtDNA OXPHOS genes achieved significance, particularly those in complex IV. In addition, several genes were consistently differentially expressed across nuclear backgrounds, including CG34092 (ND1), CG34063 (ND2), CG34076 (ND3), CG34086 (ND4L), CG34090 (CytB), and CG34073 (ATPase6).

In females the patterns were different and a much larger suite of OXPHOS-associated genes were differentially expressed from both mtDNA and nuclear-encoded genes. This effect was mainly isolated to the AutW132 nuclear background, which demonstrated a high degree of DE driven by alternative mtDNAs. In AutW132, genes that were differentially expressed in males were more significantly differentially expressed in females. The opposite effect is evident in the OregonR nuclear background, in which only a small number of OXPHOS-associated genes were differentially expressed in females.

Global gene expression and mtDNA substitution

The different patterns we observed between males and females in the OXPHOS gene subset suggest the relationships between gene expression and mtDNA substitution differs between the sexes. To investigate whether this effect is evident at a global gene-expression level, we plotted the log2-fold change in gene expression between siI and OregonR mtDNAs for each nuclear background within each sex (Figure 5). We found there were dramatic differences between the sexes in the patterns of log2-fold change; the majority of genes were significantly positively correlated in females (global gene correlation r = +0.14, t = 16.42, d.f. =12,623, P < 2.2e−16), but negatively correlated in males (r = −0.12, t = −13.99, d.f. = 13,165, P < 2.2e−16). Notably, the mtDNA genes (highlighted in Figure 5) were consistently positively correlated between nuclear backgrounds in both sexes (females: r = +0.97, t = 12.56, d.f. = 11, P = 7.279e−08; males: r = +0.86, t = 5.51, d.f. = 11, P = 0.0002). Figure 5 shows the results for the unfiltered data set. The plots of the filtered data are shown in Figure S4. In both sexes, the genes from OXPHOS complex IV and ATP synthase are generally clustered together and this may reflect their stoichiometric dependence and strict regulation. Complex I genes, by contrast, demonstrate wide variation in log2-fold change as a result of mtDNA substitution. The CG34063 (ND2) gene is the only gene that is consistently upregulated in both nuclear backgrounds and in both sexes as a result of mtDNA substitution. For the mtDNA OXPHOS genes, the ND2 gene was consistently ranked as the highest upregulated log2-fold change and ND3 was consistently the lowest ranked (highest negative log2-fold change).

Figure 5

Mitonuclear effects on gene expression differ between the sexes. Biplots of mtDNA effects on gene expression are shown for (A) females and (B) males. The effects of mtDNA substitution are reported as log2-fold changes in the OregonR and AutW132 background on the abscissa and ordinal, respectively. mtDNA genes are labeled with their gene identifiers and show positive correlations between nuclear backgrounds in both sexes. Nuclear OXPHOS genes (n = 73) are shown as red points. ND2 and ND3 genes were consistently ranked as the highest and lowest mtDNA genes, respectively. mtDNA genes are color coded by complex: complex I, black; complex III, purple; complex IV, green; ATP synthase, gray. Global gene correlations are reported in the main text.

Clustering genes by between-genotype expression profiles

Using a clustering approach, we identified subsets of genes from the global distribution that demonstrated distinct expression patterns across the four genotypes. Theoretical profiles corresponding with mtDNA, nDNA, and mtDNA × nDNA interactions are described in Figure S1. In total, we produced k = 20 clusters for each sex and the filtered data set was used, since a large number of genes in the unfiltered data have zero read counts. Hybrid tree topologies and clusters are shown for females and males in Figure S5 and Figure S6, respectively. We identified a number of clusters that were consistent with first-order nuclear effects in both females and males (Figure 6 and Figure S7, respectively). Clear mtDNA effects were only apparent in the female data set and males showed very little evidence of mtDNA effects (see above).

Figure 6

Gene clusters demonstrating a spectrum of genetic effects in females. The abscissa shows the (mito; nuclear) genotype. The log-fold change is shown on the ordinal as determined by MBClusterSeq. The black lines outline individual zero-centered gene profiles across genotypes. The red line is the per-genotype mean value across all genes in the cluster. The cluster figure for males is shown in Figure S7. Colored squares show the main genetic effect captured by the cluster, as cartooned in Figure S1: red, nuclear effect; blue, mtDNA effect; green, nuclear + mtDNA effect; black, mtDNA × nDNA interaction. A, AutW132; O, OregonR.

Males showed an overwhelming enrichment of mitonuclear interactions across gene clusters (Figure S7). The effects of mtDNA substitution are largely in opposite directions in the alternative nuclear backgrounds, an effect that confirms the patterns of global gene expression we observed in Figure 5 using an independent analysis tool. In contrast, the female data set showed only a few examples of clusters consistent with mitonuclear interactions (black squares in Figure 6).

Cluster GO-enrichment analysis

To better understand which GO processes are associated with mitonuclear interactions, we performed a GO-enrichment analysis on each cluster in each sex. In females, 16/20 clusters contained GO terms that were significant (P < 0.05) after Bonferroni correction (Antonov ), whereas only 13/20 clusters contained enriched terms in males. The full lists and significance of the GO terms are shown in Figure S8 (females) and Figure S9 (males).

Clusters with clear mitonuclear effects in females were clusters 17 and 18 (Figure 6). The most significantly enriched GO terms in these clusters were generally related to mitochondrial-ribosomal and protein translational processes (see clusters 17 and 18 in Figure S8).

edgeR statistical analyses

The clustering approach provided qualitative support for an excess of mitonuclear interactions in males, and a larger proportion of mtDNA effects in females. To formally test the effects of mtDNA, nuclear, and mitonuclear variation, we conducted a DE analysis using edgeR (Robinson ) in each sex separately. We conducted separate tests in males and females because the RNA-seq libraries have inherent differences in dispersion (biological coefficients of variation) between sexes. Stratifying the analyses by sex ensured that we made comparisons of genetic effects using appropriate dispersion estimates without conflating the DE estimates.

Figure 7, A and B, describes the distributions of significantly differentially expressed genes (P < 0.05) in a sex-biased expression context. The plotted data are the mean expression values across all RNA-seq libraries in each sex. Significantly, DE genes corresponding to the three forms of genetic variation (mtDNA, nuclear, and mitonuclear) are highlighted as black, purple, and green data, respectively. In females (Figure 7A), there was a large overlap between DE genes from all three sources of genetic manipulation but these significant genes were not enriched in female- or male-biased genes. Likewise in males, the majority of genes significantly differentially expressed by nuclear type were found in both non-sex-biased regions of the distribution, along with male-limited regions that correspond with testes and sperm-related biological processes (Figure 1, Table 1). While we found a few genes in the male-limited region of the distribution that were DE by mtDNA variation, DE genes associated with mtDNA haplotype variation were not enriched in this region (Figure 7B).

Figure 7

The distribution of significant mtDNA, nDNA, and mitonuclear genotypes in expression space. Results of edgeR analyses on the complete data set are shown. Female expression is plotted on the abscissa, males on the ordinal. The position of differentially expressed genes by nuclear (black), mtDNA (purple), and mtDNA × nDNA (green) variation are shown for (A) females and (B) males. The absolute numbers of significant genes in each category are described in bar plots in (C). Results for each category are divided into up- and downregulated genes. For each category females are shown leftmost, males rightmost.

Across sexes, there were magnitudinal differences in the numbers of genes differentially expressed by mtDNA, nuclear, and mtDNA × nuclear variation (Figure 7C); females showed larger responses to mtDNA and mitonuclear interactions, whereas males demonstrated marginally greater numbers of nuclear DE genes. For mtDNA variation (siI/OregonR), females demonstrated upregulation of 70 genes, and 39 genes were upregulated in males. There were 621 downregulated genes in females and 26 in males. By far the largest source of DE was for nuclear genetic variation (AutW132/OregonR) in both females (1737 up, 1952 down) and males (1836 up, 2246 down). Far fewer genes were associated with mitonuclear interactions (contrast siI; AutW132 vs. OregonR; OregonR) and these included 1462 in females and 556 in males. Interestingly, males showed larger numbers of DE genes in the mitonuclear category than in the first-order mtDNA category, consistent with the patterns we observed in the clustering analysis (see above).

To determine if genes behaved consistently across the sexes in response to genetic manipulation, we intersected those genes that were significantly differentially expressed (P < 0.05) by mtDNA, nuclear, and mitonuclear variation from females and males. The between-sex intersections are shown in Figure 8. There were 26 genes that were consistently differentially expressed by mtDNA variation, 1426 genes differentially expressed by nDNA variation, and 124 that were consistently differentially expressed by mitonuclear variation (Figure 8; gene identifications of the intersected genes can be found in Figure S10, Figure S11, and Figure S12). mtDNA genes are enriched in the intersected genes for mtDNA and nuclear variation, and ATPase6, ND5, and ND6 demonstrate a consistent mitonuclear effect across the sexes. We found a large overlap between the consistently differentially expressed genes identified by DESeq (Table 2) and those identified by edgeR (Figure S10) as a consequence of mtDNA variation.

Figure 8

Robust differentially expressed genes across (A) mtDNA, (B) nDNA, and (C) mtDNA × nDNA categories. Venn diagrams describe the intersection of gene identifiers by sex in each category, as determined by DE analysis in edgeR on the complete data set. The area of the circles is relative to the number of genes within a category. Females are shown on the left, and males on the right of each diagram. Genes within the intersection between females and males are listed in Figure S10, Figure S11, and Figure S12 for each category, respectively.

Validating RNA-seq using qPCR

We used six nuclear genes that showed DE by nuclear variation and one internal control gene (rp49) to test whether the fold changes we measured using RNA-seq were correlated with qPCR. The methods and results can be found in File S1. The log2-fold change between AutW132 and Oregon R in each mtDNA background was calculated for each gene in both sexes using the comparative threshold cycle (CT) method (2−ΔΔCT) (Schmittgen and Livak 2008). We found a significant positive correlation between the log2-fold changes measured with RNA-seq and qPCR (r = +0.69, t = 4.52, d.f. = 22, P = 0.00016; Figure S13).

Discussion

We tested for the effects of mtDNA, nuclear, and mitonuclear variation on gene expression variation and found both genotype- and sex-specific responses to all three forms of genetic variation. Females showed overall stronger effects of genetic variation than males, and males showed little evidence of sex-limited expression of DE genes associated with mtDNA variation. mtDNA variation was most significantly associated with variation in mitochondrial OXPHOS gene expression and particularly those genes that are encoded by the mtDNA. In this panel of genotypes we found an overall opposing direction of mtDNA variation on global gene expression between the sexes, yet a conserved effect of mtDNA-encoded genes. Our findings are important to the Drosophila community, in which our reference strain, OregonR, has been widely adopted as a commonly used wild-type laboratory stock. Ongoing work in our laboratory aims to understand whether these effects are universal across multiple nuclear backgrounds, and not restricted to this particular nuclear OregonR-AutW132 pair. We discuss our results in the context of OXPHOS protein organization, mtDNA–nDNA covariation, and the evolutionary significance of sex and genotype interactions.

Mitochondrial gene expression and fitness

Gene expression variation in mtDNA and nuclear-encoded mitochondrial genes has been associated with a large number of human diseases, including cancer (Penta ), neurodegenerative diseases (Schapira 1998; Swerdlow and Khan 2009), aging (Tońska ), and type 2 diabetes mellitus (Mootha ); however the exact mechanisms of action are poorly understood (Borowski ). Early studies suggested the expression of OXPHOS genes in early development conditions the rate of electron transport enzyme activity throughout life in Caenorhabditis elegans (Dillin ). More recently, knockdown of mitochondrial ribosomal proteins causes mitonuclear protein imbalance, reduced respiration, and activation of the mitochondrial UPR (Houtkooper ). Paradoxically, this imbalance confers an overall positive, hormesis-like effect on life span and appears to be conserved between C. elegans and mammalian (mouse hepatocyte) cell lines. The regulation of mitochondrial genes is therefore important for organismal health, and transcript (or protein) imbalance between mitochondrial proteins may provide one arena for a cell to sense aberrant gene or protein function.

mtDNA substitution affects mtDNA expression

We found expression levels of mtDNA genes to be among the most differentially expressed between mtDNA haplotypes. Importantly, these differences were not unidirectional across all genes. In all nuclear background × sex combinations, the siI haplotype conferred upregulation of expression in some genes (e.g., ND2 gene), whereas genes from the same complex (complex I, e.g., ND3) were among the most downregulated genes as a consequence of mtDNA substitution. This suggests that the relative abundance of transcripts is lower in siI haplotypes, however, there are exceptions to this rule (e.g., ND2). For mtDNA genes, we found a strong positive correlation between log2-fold changes in each nuclear background, suggesting that mtDNA genes behave similarly in different nuclear backgrounds and sexes. Importantly, there was a consistent rank order for the genes that were most up- and downregulated as a consequence of mtDNA substitution.

Mitochondrial haplotype variation has previously been shown to modify mtDNA copy number variation and mtDNA protein coding gene expression on a common nuclear background in Drosophila (Camus ). In that study, the ND5 gene showed expression values consistent with phenotypic divergence, suggesting the mtDNA polymorphisms may behave as expression quantitative trait loci. The ND2 gene was not assayed in Camus et al.’s study due to its high A+T content, so comparisons between the present study and that study cannot be made. We show that mtDNA substitution can have a large impact on relative gene expression when large numbers of nucleotide polymorphisms are present in the contrasting genotypes (e.g., ∼103 amino acid substitutions and 418 synonymous SNPs between Oregon R and siI haplotypes). The results from the current study show that the impact of mtDNA haplotypes on mtDNA gene expression may be an additive effect of the point mutations that differ between the compared haplotypes.

Protein sequences or regulatory sequences?

Using the mitonuclear introgression model we specifically substituted alleles in the mtDNA coding region. We do not know the extent of mutations in the regulatory control region of the mtDNAs (D-loop), which is where transcription is initiated in D. melanogaster (Goddard and Wolstenholme 1978) and D. simulans (Goddard and Wolstenholme 1980). In both species, at least one (Clayton 1982) origin of replication occurs roughly in the center of the A+T-rich region (Goddard and Wolstenholme 1978; Wolstenholme 1992; Lewis ; Torres ). Our findings suggest that protein products of genes that assemble in a single complex (complex I) can have drastically different transcript abundance, and the mutations that segregate between siI and OregonR haplotypes alter expression in a gene-by-gene basis.

Complex I is one of the largest and most complicated enzymes in the eukaryotic cell (Vinothkumar ), and its crystal structure has been resolved across a wide range of species, e.g., Escherichia coli (Efremov and Sazanov 2011), Thermus thermophilus (Efremov ), and Bos taurus (Vinothkumar ). Among the consistent characteristics are the protein identifiers and their spatial organization within the membrane domain (Vinothkumar ). One of the earliest genome-wide coexpression analyses provided good evidence in Saccharomyces cerevisiae that genes that physically interact or are present in the same metabolic pathway have similar expression levels (DeRisi ). Our results show that a “healthy” fly can exhibit wide variation in mtDNA transcripts from some complexes (e.g., complex I), yet consistent coexpression patterns across other complexes (e.g., complex IV and ATP synthase). Taken together, these results suggest some transcripts are more sensitive to genetic polymorphism than others, and that transcript variation within a protein complex may provide a useful trait to characterize the effects of mutations, and whether these effects are associated with deleterious phenotypes. Finally, the phenotypic effects on separate mitochondrial respiratory functions could help pinpoint whether complexes with tightly regulated transcripts perform better than those that have wide variation. If wide variation is deleterious, we would predict that complex I function, for example, would show more wide-ranging phenotypes between siI and OregonR haplotypes, than complex IV.

Why would ND2 be massively upregulated and ND3 massively downregulated as a result of siI mutations? In spite of their physical proximity in the mitochondrial inner membrane domain, we found those OXPHOS genes that were most sensitive to mtDNA variation are adjacent proteins (Vinothkumar ). For example, Figure 5 shows that gene transcripts ND3 and ND4L are closely grouped in the plot. Likewise, ND2 and ND4 are consistently upregulated in both nuclear backgrounds and in both sexes. The patterns of expression are coextensive with the position of protein products; ND3 and ND4L, and ND2 and ND4 are adjacent proteins in the membrane. Given that there are only seven mtDNA-encoded proteins in the membrane domain it is possible that this observation would be expected by chance, however it does suggest that conversion rates of transcript to protein may depend on the position of a protein in the membrane.

We have previously identified a number of SNPs that are different between the coding regions of siI and OregonR haplotypes and which have putative deleterious effects on protein function (Mossman ). One of the private mutations to the siI-OregonR pair resides in the ND2 gene and is two amino acids downstream of a putative deleterious amino acid polymorphism. Future work will aim to characterize whether this putative mutation is involved in the upregulation of the gene, possibly as a response to deleterious protein function; resulting in overall genetic robustness. Previous studies have shown that gene transcription responds to genetic mutation and may rescue phenotypes via a compensatory network (Kafri ; Rossi ). Moreover, during mammalian mtDNA transcription, polypeptides are theoretically transcribed in a 1:1 ratio because the polycistronic transcripts are fully transcribed from both mtDNA duplex strands (Gagliardi ). In Drosophila mtDNA there are five different transcription initiation sites and transcript cleavage is signaled by the presence of cloverleaf structures of tRNAs (Ojala ). There is considerable heterogeneity in mRNA expression between mtDNA genes within OXPHOS complexes, and these differences are largest in complex I (NADH dehydrogenase) (Torres ). We also found NADH dehydrogenase to have the widest variation in log2-fold change when comparisons were conducted within a gene (as a consequence of mtDNA variation). That is, mtDNA genetic variation was associated with differences in log2-fold change. Taken together, these results suggest complex I gene expression shows wide variation (i) across genes within a complex and (ii) within genes as a consequence of mtDNA variation.

Following cleavage via endonucleases, polypeptide mRNAs are processed with a theoretical uniform abundance, which we did not observe. Given our haplotypes were from separate species, and those species have different control regions, with variable length and possibly SNP polymorphisms, it is possible that these regions alone are associated with the expression differences, with SNP variation having negligible effect. However, we favor the hypothesis that SNP variation does impact transcript levels, simply because a main effect of the regulatory (D-loop) region would presumably affect transcript levels in a more universal way, with little evidence of gene-specific expression patterns. Heterogeneity in transcript abundance across genes in Drosophila has been suggested to result from various post-transcriptional mechanisms including differential transcript stability and differences in the processing of mature transcripts (Torres ). Our results suggest mtDNA variation can modify mtDNA gene transcriptional regulation and that complex I is particularly sensitive to mtDNA variation in this mtDNA haplotype pair (OregonR vs. siI).

Gene-by-gene interactions for transcription

mtDNA-encoded protein genes showed the most consistent patterns of DE identified using multiple analysis tools. Three mtDNA genes (ND5, ND6, and ATPase6; Figure S12) demonstrated a significant mitochondrial × nuclear effect in both sexes. That is, the effect of mtDNA substitution in one nuclear background was different from the effect in the other. As expected, nDNA variation was associated with DE in many genes in both males and females, however, we found females were generally more sensitive to mtDNA variation, particularly for OXPHOS complex genes. While males showed gene clusters consistent with mitonuclear effects, these were generally low magnitude differences and were not detected using DE analyses. We also found a larger effect in the AutW132 nuclear background than in the OregonR nuclear background in females. Our gene lists (in Figure S10, Figure S11, and Figure S12) documenting the genes that were consistently differentially expressed across sexes are likely conservative estimates. On the other hand, they represent core sets of genes that are robustly differentially expressed across genetic backgrounds and sexes. Our clustering analyses showed clear evidence that many transcripts have gene-by-gene interaction patterns in males, yet mainly mtDNA effects in females; a pattern that was consistent across analysis tools. Interestingly, the nuclear genes involved in OXPHOS showed mitonuclear patterns in males, but no mitonuclear effects in females (red dots in Figure 5). Using an agnostic clustering approach, we identified sets of genes that showed similar expression patterns across genotypes and which were enriched for GO categories associated with mitochondrial function. It is now well established that genes in the same pathway share similar levels of expression (Stuart ; Kafri ) and we detected a strong signal that mitonuclear interaction gene clusters were enriched for processes involving translation and mitochondrial metabolism in both sexes. Moreover, the mtDNA OXPHOS genes tended to cluster in complex-specific clusters and this effect was more noticeable in females than in males. Similar results have been observed in humans where genes from distinct OXPHOS complexes tend to cluster together (van Waveren and Moraes 2008).

The female transcriptome is more sensitive to genetic variation

Considering the global gene set, we found that females exhibit clear genotype-specific transcript profiles regardless of the gene set analyzed. In contrast, the MDS plots revealed males are distinguishable by nuclear backgrounds across all genes, but the effects of mtDNA are only present in the most differentially expressed gene sets and mainly restricted to the AutW132 nuclear background. This first suggests that mtDNA genes are among the most differentially expressed in males, and second, that mtDNA × nDNA interactions are present in males, and these operate on a large subset of the complete gene set. These were generally smaller in magnitude than in females. This phenomenon may help explain why a greater number of mitonuclear interactions were observed in males via gene clustering and these made detection of first-order effects of mtDNA and nuclear variation more difficult as fewer genes achieved statistical significance with DE-analysis tools. Using a combination of DE analyses provides clearer evidence that male and female transcriptomes respond differently to mitonuclear variation, and suggests a much larger number of genes are sensitive to mitonuclear variation in males, even though these are not “significantly” differentially expressed.

Retrograde signaling is more prevalent in females

We found mtDNA substitution in females conferred large effects on nuclear genes, suggesting pervasive retrograde signaling (mtDNA-to-nDNA feedback) between the genomes, and a clear separation of genotypes based on transcript expression. The effects of mtDNA variation were less obvious in males and a significantly smaller number of nuclear genes responded to mtDNA variation. While we acknowledge that there could be inherent differences between the male and female RNA-seq libraries due to a block effect, we conducted the majority of analyses within a sex to minimize a block effect on analyses interpretation. The magnitude and direction of log2-fold change in mtDNA genes were similar across sexes, suggesting we were able to capture real biological signal for the most differentially expressed genes. Importantly, when correlations (Figure 5) were compared across sexes, our mitonuclear results demonstrate a clear sexual dimorphism, even though the within-sex effects were calculated with sex-specific dispersion parameters. Under a null hypothesis of no sex-specific mitonuclear effect, we would not expect the sexes to differ in the direction of effect for the global pattern of expression. Given that the majority of male transcripts did not respond significantly to first-order mtDNA variation, the global transcript negative correlation we observed between nuclear backgrounds may be a product of reduced signal of mtDNA effects. However, mtDNA genes did respond and the mtDNA effect, while small for the majority of the transcriptome, was in the opposite direction in females. Our results show that the sensitivity to detect first-order effects of mtDNA variation are hindered by mitonuclear interactions, and a clustering approach is more likely to capture interaction effects. One interesting question arises: is this evidence of sexual antagonism (Rice 1984) for gene expression?

Genes from both mtDNA and nDNA genomes can respond to selection in females. The female environment likely provides the only opportunity for mitonuclear expression to be exposed to selection, even though genes spend, on average, half their lifetime in each sex. Under this scenario, it has been suggested that nuclear genes that interact with mtDNAs would benefit from being on the X chromosome because there is greater opportunity for cotransmission and therefore coadaptation with mtDNA (Rand ; Wade and Goodnight 2006). In Drosophila, the X chromosome has been implicated in mitonuclear epistatic interactions (Rand ; Montooth ), providing an opportunity to reduce intraindividual genetic conflict (sensu Werren 2011). However, mitonuclear genes are underrepresented on the X chromosome in Drosophila (Rogell ), a pattern that is consistent also across various mammals (Drown ) and C. elegans (Dean ). Further, nuclear-mitochondrial gene duplications, which could help mitigate sexual conflict, are rarely relocated to the X chromosome in Drosophila (Gallach ). In contrast, birds do not show under- or overrepresentation of mitonuclear genes on the Z chromosome (in birds, males are the homogametic sex and genes on the Z chromosome were tested; Drown ).

Alternative sexual conflict-based hypotheses have been proposed to explain these observations (Drown ; Dean ; Rogell ). Our results suggest pervasive retrograde signaling occurs between mtDNA and nuclear genes, motivating that future investigations on mitonuclear coevolution should focus on key regions of the nuclear genome that harbor mitonuclear interacting genes, and not just canonical OXPHOS or mitochondrial genes encoded by nDNA. For example, do these nuclear loci demonstrate evidence of unique patterns of selection compared to closely linked loci? We identified two nuclear genes, and Jonah 25Bi, that show consistent evidence of retrograde signaling in all genotypes and in both sexes. These genes are sensitive to mtDNA polymorphism, yet are not OXPHOS genes. We have also described a gene list of mitonuclear genes that are attractive targets of such a study (Figure S12).

mtDNA effects are not limited to males or male-limited genes

The Frank and Hurst hypothesis (Frank and Hurst 1996)—that males are more sensitive to mtDNA substitution—has limited support in this study. We found the mtDNA protein coding genes to be enriched in the gene list that was differentially expressed by mtDNA variation, and this effect was consistent across nuclear backgrounds. We also found a similar effect in females. While these findings suggest males are not suffering more than females as a result of mtDNA substitution, the enrichment of mtDNA genes per se is some evidence that males are sensitive to exclusively OXPHOS genes (Figure 4). It is possible that the nucleotide substitutions between siI and OregonR do not manifest with sufficiently large effect sizes on globally expressed genes to be detected in males, as has been observed in other mtDNA-variation studies on gene expression (Innocenti ). Moreover in Innocenti et al.’s study, the effects of mtDNA substitution in a w1118 D. melanogaster background were enriched in male-limited genes. Their finding evidenced a sex-specific selective sieve, whereby mtDNA mutations that were not under selection in males (due to maternal inheritance) could principally manifest in male-limited tissue as males are at a genetic dead end for mtDNA evolution. To independently test for mtDNA-associated expression differences in sex-limited genes, we also conducted our investigation on whole flies, including all reproductive tissue. We acknowledge that whole fly analyses cannot resolve sex-biased expression that may be due to differences in tissue contributions. However, we wanted to make our analyses comparable to other whole organism gene expression studies in Drosophila. We found no differences in fecundity between the strains used in this study (see File S1, Figure S14, Table S2, Table S3), suggesting the rate of turnover of gametes, and hence reproductive tissue mass, was not significantly different between strains and was unlikely to bias gene expression values between genotypes. Pathological mtDNA mutations principally manifest in tissues with high ATP demand (reviewed in Stewart and Chinnery 2015). The tissue-specific pathological effects of mtDNA variants in humans is a good motivation for future work to specifically test gene expression differences in tissues with high ATP demand in fruit flies, such as neuronal tissue, or flight muscle.

We did not recapitulate the sex-specific selective sieve effect in our study, in spite of (i) significant numbers of genes that demonstrated sex-specific gene expression (Figure 1), and (ii) large amounts of nucleotide variation between our haplotypes (>100 amino acid substitutions, and >400 synonymous mtDNA polymorphisms). Instead, we found no enrichment for mtDNA effects in male-biased genes [in regions of the gene expression space associated with male-limited traits (e.g., spermatogenesis and testes)]. We only considered a pairwise comparison of mtDNA variation, in contrast to Innocenti et al.’s five haplotypes, and hence do not have the same polymorphisms present in our experiment. However, it is possible that the w1118 nuclear background used in Innocenti ) is more sensitive to mtDNA variation than other nuclear backgrounds and these differences may manifest disproportionately in males. In the present study we found greater mtDNA sensitivity in the AutW132 background, providing good evidence that nuclear background per se can alter the sensitivity of the transcriptome to mtDNA variation. A main take home message in the present study is that there are differences between nuclear backgrounds and these can modify the effects of mtDNA variation. Extensive evidence of mtDNA × nDNA interactions on phenotypes (Mossman ), including comparisons with the w1118 background (Zhu ), suggests mtDNA haplotype-associated phenotypes are not always general results across nuclear backgrounds. We are now assessing the effects of environment on mitonuclear interactions to test whether mitonuclear effects are influenced by environment (gene-by-gene-by-environment interactions), or whether genotypes are robust to environmental perturbation. Previous studies suggest environment can have a large and sometimes unpredictable effect on mitonuclear interactions (Mossman ), however, there are no studies quantifying this in Drosophila gene expression.

In summary, we found considerable variation in transcript expression as a result of mtDNA, nDNA, and mitonuclear substitution. Contrary to the Frank and Hurst hypothesis, we did not find enrichment of mtDNA effects targeting male-limited gene expression in males. However, we did find large sex differences in the effects of genetic manipulation. In general, females showed larger genetic effects on transcript abundance and female genotypes were distinguishable based on their global gene-expression patterns, inconsistent with the lower sensitivity in males. In males, the majority of transcripts demonstrated mitonuclear effects in clustering analyses. Although mtDNA genes were preferentially differentially expressed, we interpret this as evidence that both mtDNA and nDNA covariation are important for transcript expression in this genotype panel. Future work will identify: (i) if these effects are general over a larger suite of genotypes, and (ii) by what mechanisms of action do mtDNA haplotypes affect gene expression. For example, are the polymorphisms in genic regions of the mtDNA responsible, or is variation in the mtDNA regulatory region key to the observed expression patterns? Are these effects tissue specific? Our general findings suggest that mtDNA effects can vary between nuclear genetic backgrounds depending on the sex tested, and therefore therapeutic methods to overcome mitochondrial diseases in humans should consider mitonuclear covariation as potential sources of phenotypic variation and therapy outcomes.