Mitonuclear mismatch alters nuclear gene expression in naturally introgressed Rhinolophus bats

Abstract

Background

Mitochondrial function involves the interplay between mitochondrial and nuclear genomes. Such mitonuclear interactions can be disrupted by the introgression of mitochondrial DNA between taxa or divergent populations. Previous studies of several model systems (e.g. Drosophila) indicate that the disruption of mitonuclear interactions, termed mitonuclear mismatch, can alter nuclear gene expression, yet few studies have focused on natural populations.

Results

Here we study a naturally introgressed population in the secondary contact zone of two subspecies of the intermediate horseshoe bat (Rhinolophus affinis), in which individuals possess either mitonuclear matched or mismatched genotypes. We generated transcriptome data for six tissue types from five mitonuclear matched and five mismatched individuals. Our results revealed strong tissue-specific effects of mitonuclear mismatch on nuclear gene expression with the largest effect seen in pectoral muscle. Moreover, consistent with the hypothesis that genes associated with the response to oxidative stress may be upregulated in mitonuclear mismatched individuals, we identified several such gene candidates, including DNASE1L3, GPx3 and HSPB6 in muscle, and ISG15 and IFI6 in heart.

Conclusion

Our study reveals how mitonuclear mismatch arising from introgression in natural populations is likely to have fitness consequences. Underlying the processes that maintain mitonuclear discordance is a step forward to understand the role of mitonuclear interactions in population divergence and speciation.

Background

Mitochondria are crucial to the fitness of an organism by regulating energy production in the cell. In particular, the oxidative phosphorylation (OXPHOS) pathway, necessitates epistatic interactions between proteins encoded by both mitochondrial and nuclear (mitonuclear) genomes [1, 2]. Because mitochondrial DNA (mtDNA) typically has a higher mutation rate than nuclear DNA [3], those nuclear genes that code for mitochondrial proteins (N-mt genes) evolve compensatory mutations in order to maintain compatibility of mitonuclear genomes [4, 5]. Selection for this mitonuclear compatibility leads to co-evolution of these two genomes, resulting in mitonuclear co-adaptation in different isolated populations [6,7,8].

The co-adaptation between mitochondrial and nuclear genomes can be disrupted in cases where formerly isolated populations come into secondary contact and undergo hybridization with introgression of mtDNA [9]. In particular, in cases of extensive introgression of mtDNA with little or no introgression of nuclear DNA (nDNA), the mitochondrial genome becomes mismatched with its nuclear genetic background [10], potentially resulting in reduced fitness due to mitonuclear conflicts and/or incompatibility [11,12,13,14,15]. Therefore, mitonuclear conflict caused by mitochondrial and nuclear DNA mismatch (mitonuclear mismatch) has been documented as an important driver of intrinsic reproductive isolation [16] and may play an important role in speciation [17, 18].

Mitonuclear mismatch has been observed in many naturally occurring populations [19,20,21] and also in humans [22], and is commonly detectable as mitonuclear discordance [23]. Studying natural individuals that have identical or similar nuclear genetic backgrounds but with markedly different mitochondrial variants can help to understand the nature of mtDNA-nDNA interactions and the evolutionary consequences of mitonuclear mismatch on population divergence. So far, few studies have directly investigated the impact of mitonuclear mismatch on fitness in natural populations (but see [10, 24]), largely because it is difficult to compare the fitness of individuals with different combinations of mitochondrial and nuclear genotypes. To tackle this limitation, we can use indirect ways of assessing possible fitness consequences, such as by examining transcriptional differences in multiple tissues between individuals with matched and mismatched mitonuclear genomes. Up to now, only a few studies have used comparative transcriptomics to investigate the effects of mitonuclear interactions on gene expression, and nearly all such studies focused on the model organism Drosophila [25,26,27]. Most of these studies indicated that mtDNA variation could alter nuclear gene expression, particularly when individuals with mismatched mitonuclear genomes are placed under stress (e.g. hypoxic conditions). As far as we know, few such studies have been performed in naturally introgressed populations [28, 29], which are critical to understand how individuals with mismatched mitonuclear genotypes have evolved to cope with mitochondrial dysfunction, such as reduced ATP production [30] and elevated oxidative damage [31, 32].

Bats are the only mammals to have evolved powered flight, and, for this reason, they have extremely high metabolic rates compared to other similar-sized mammals [33]. It is therefore likely that bats are unable to tolerate mitochondrial dysfunction caused by mitonuclear mismatches. Here we focus on the intermediate horseshoe bat (Rhinolophus affinis) which includes three recently diverged subspecies in China: two mainland subspecies (R. a. himalayanus and R. a. macrurus) and one island species from Hainan (R. a. hainanus) (Fig. 1a, see also [34, 35]). These three subspecies diverged less than one million years ago, and our previous phylogenetic studies of this species have indicated that himalayanus colonized Hainan to form hainanus during a period of glaciation, which then recolonized the mainland to form macrurus [34,35,36]. Extensive mitochondrial introgression (almost complete replacement of mitochondrial genome) with little or no nuclear introgression, appears to have occurred from macrurus to himalayanus in regions of secondary contact, following the recent population expansion of macrurus soon after the last glacial maximum (~ 18,000 years ago) (Fig. 1a, b, see also [34,35,36,37]). The discordant mitonuclear introgression between these two subspecies has resulted in mismatched mitonuclear genomes of R. a. himalayanus individuals, providing a unique model system in which to investigate the transcriptional impact of mitonuclear mismatch in naturally introgressed populations. Here we compare individuals with matched mitonuclear genotypes (both nuclear genome and mitogenome are from R. a. himalayanus) to ones with mismatched mitonuclear genotypes (nuclear genome is from R. a. himalayanus while mitogenome is from R. a. macrurus) (Fig. 1c).

Fig. 1

Experimental design to study the transcriptional effects of mitonuclear mismatch in R. affinis. a Phylogenetic relationships between the two hybridizing subspecies, R. a. himalayanus and R. a. macrurus. Extensive mtDNA introgression has occurred from macrurus to himalayanus in their secondary contact regions. b Geographic distribution of macrurus and himalayanus and the sampling locality of himalayanus in this study. c The study system includes five males of each himalayanus group, mitonuclear matched (both nuclear and mitochondrial genomes are from himalayanus) and mismatched (nuclear and mitochondrial genomes are from himalayanus and macrurus, respectively) groups. d We sampled six adult tissues (pectoral muscle, heart, brain, liver, cochlea, and small intestine) to assess transcriptional variations of nuclear genes between matched and mismatched groups

Because mitochondria in mitonuclear mismatched individuals tend to generate higher levels of reactive oxygen species (ROS) [38, 39] and elevated oxidative damage [31, 32], we hypothesized that genes with roles in the protection of cells against oxidative stress would be significantly upregulated in mismatched individuals. To test this hypothesis, we assess transcriptional variations of nuclear genes in multiple adult tissues (Fig. 1d). Specifically, we examine differentially expressed genes (DEGs) between mitonuclear mismatched and matched individuals (also see [25,26,27]). These DEGs may be essential for cells or organisms to respond to the cellular stress caused by mitonuclear mismatch in natural conditions.

Methods

Sampling and experimental design

We previously found that many Rhinolophus affinis himalayanus individuals sampled from Anhui and Jiangxi provinces showed introgression of mtDNA from R. a. macrurus [37]. In this study we captured 10 adult males of R. a. himalayanus from one cave in Anhui province, China in 2019 (Fig. 1b). R. a. himalayanus is not an endangered species in the local region. Bats were euthanized by cervical dislocation. For each individual, we collected six tissues (pectoral muscle, heart, brain, liver, cochlea, and small intestine) with RNase-free tubes. We also collected the brain tissue from four R. a. macrurus individuals used in our previous study [40] in order to examine the sequence differences between the two subspecies (see below). Tissues were frozen immediately in liquid nitrogen and stored in a − 80 °C freezer.

To identify himalayanus individuals with or without introgression of macrurus mtDNA, we reconstructed phylogenetic relationships between these newly sampled individuals and previous ones from macrurus and himalayanus. For this, we sequenced the cytochrome b (cytb) gene from new samples and downloaded cytb sequences published previously from five macrurus and five himalayanus [36]. Primers and PCR profiles for amplification of cytb gene have been described previously [41]. MEGA 6 [42] was used to align all sequences and a Neighbor-Joining (NJ) tree was reconstructed with 1000 bootstraps. The NJ tree revealed that five of the new himalayanus individuals had mitochondrial haplotypes that clustered with macrurus while the remaining five individuals clustered with himalayanus (Additional file 1: Supporting Figure S1a). Thus, we separated these 10 new himalayanus individuals into two groups, one of five individuals with macrurus mtDNA (mismatched group, Nc-himalayanus:Mt-macrurus) and the other one of five individuals with himalayanus mtDNA (matched group, Nc-himalayanus:Mt-himalayanus) (Fig. 1c and Table 1).

Table 1 Details of 10 R. a. himalayanus individuals used for RNA-seq analyses. Among them, two individuals (Him-16 and Him-17) used for whole-genome resequencing analysis were shown with *

Transcriptome sequencing, genome resequencing and raw data trimming

For each tissue sample, we extracted total RNA extraction with TRIzol (Life Technologies Corp., Carlsbad, CA, USA) and conducted library constructions using Illunima’s TruSeq mRNA Standard library preparation kit. One library from the muscle tissue was discarded due to its low quality. All the remaining 63 libraries were quantified and qualified with Agilent 2100 Bioanalyzer and sequenced with Illumina HiSeq X Ten (paired-end 150 bp). Overall, we generated RNA-seq data for 59 himalayanus and four macrurus samples and obtained an average of 22 million reads per sample (Additional file 4: Table S1).

We generated ~ 30 Gb data with the whole-genome sequencing coverage of 15 × for each of the two individuals (Him-16 and Him-17) representing mismatched and matched groups, respectively (see below). Genomic DNA was extracted using muscle tissue with DNeasy kits (Qiagen) and was quantified using a Qubit 2.0 Fluorometer. DNA libraries were constructed using the NEB Next Ultra DNA Library Prep kit (New England Biosciences) and sequenced on the Illumina HiSeq X Ten (paired-end 150 bp).

Raw reads from each sample were trimmed with TRIMMOMATIC version 0.38 [43] using a sliding window of 4 bp with minimum average PHRED quality score of 20 and minimum reads length of 50 bp.

Analysis of nuclear genome sequence

To characterize the nuclear genome difference between the two himalayanus groups, we first conducted an individual-based principal component analysis (PCA) with PLINK version-1.07 [44] based on SNPs datasets. For SNPs calling, we used RNA-seq reads of the brain tissue from all 10 himalayanus individuals and four macrurus. Briefly, we mapped filtered reads from each individual to the chromosome-level reference genome of R. a. himalayanus (G. Li, unpublished data) using BWA-MEM with default parameters [45]. SAMtools Version 1.9 [46] was used to generate sorted BAM file and to remove potential PCR duplicates. We used BCFtools Version 1.9 [47] to perform multisample SNPs calling. For each sample, we only kept high quality SNPs with the mapping quality (MQ) > 20, base quality (phred-scaled confidence) > 30 and sequence coverage (read depth) > 10. Second, based on SNP datasets we calculated the Weir and Cockerham estimator of Fst [48] between the two himalayanus groups using Vcftools v0.1.16 [49]. Because our main aim here is to estimate the general level of genetic differentiation between the two himalayanus groups, we applied a nonoverlapping window of 100 kb, and only windows with the minimum of 10 SNPs were used to calculate the weighted Fst values. Third, to further test for the genome similarity of the two himalayanus groups, we examined sequence differences at nuclear-encoded OXPHOS genes because their proteins show direct interactions with mitochondrial encoded proteins. We generated sequences of 68 OXPHOS genes [50]from each of the 14 individuals (10 himalayanus and 4 macrurus) (see details in Additional file 3: Supporting information and Additional file 9: Table S6). These sequences were then aligned using MEGA 6 and amino acid changes were calculated across the 14 individuals.

Analysis of mitochondrial genome sequence

We previously generated a complete mitogenome for Him-17 (GenBank accession: MT845219) [51], belonging to matched group. To compare sequence differences of the complete mitogenome between mismatched and matched groups, we generated the complete mitogenome for Him-16 (mismatched group) using similar procedures as in Ding et al. [51] (see details in Additional file 3: Supporting information). Nucleotide differences between Him-16 and Him-17 were calculated for 37 mitochondrial genes and control region using MEGA.

To investigate whether sequence differences between Him-16 and Him-17 were fixed in each group, we generated the 13 mitochondrial protein-coding genes (PCGs) for the remaining eight himalayanus individuals and four macrurus individuals (see details in Additional file 3: Supporting information). All sequences were aligned using MEGA and amino acid changes across all 14 individuals were calculated for each of the 13 PCGs. Aligned sequences of concatenated 13 PCGs were used as the input to build a Maximum-likelihood (ML) tree using RaxML with GTRGAMMA model and bootstraps 1000 setting. R. ferrumequinum (GenBank accession: KT779432) was used as the outgroup.

To complement the above nuclear analyses, we also conducted PCA and sliding-window analyses based on SNPs dataset of mtDNA using similar procedures as above. SNPs were generated using the mitogenome of Him-17 as the reference. For the sliding-window Fst analysis, we used a nonoverlapping window size of 1 Kb.

Differential expression analysis

For each tissue, filtered reads from each sample were mapped to the reference genome of R. a. himalayanus using Hisat2 [52] with default settings. Mapped reads were quantified using Featurecount [53]. We removed those lowly expressed transcripts with less than one CPM (counts per million) across samples. Read count matrices across samples were normalized in DESeq2 [54]. Prior to differential expression analysis (DE), we checked for potential outliers in each tissue using PcaGrid method [55] implemented in the rrcov R package with default parameters. Out of the 59 samples, one brain sample, two liver samples, one cochlea sample, and one intestine sample were identified as significant outliers and therefore were not used in DE analysis. Samples were renormalized after removal of outliers in each tissue. We identified differentially expressed genes (DEGs) between the two R. a. himalayanus groups in each tissue using DESeq2. A gene was considered as DEG using |log₂ (fold change)|> 1 and P value < 0.05 after Benjamini and Hochberg adjustment for multiple tests (padj < 0.05) [56].

To examine the functional role of DEGs identified in each tissue, we performed enrichment analysis using Metascape (http://metascape.org) [57]. Significance was determined using Bonferroni corrected p value of < 0.05. We reduced redundancy of significant GO (Gene Ontology) terms using the REVIGO clustering algorithm (http://revigo.irb.hr/) [58] and the clustered GO terms were visualized using scatterplots with the semantic similarities of GO terms.

Results

Mitogenome analysis

We found a total of 504 nucleotide differences (2.98% of mitogenome) between the mitogenomes of Him-16 and Him-17, which represent mismatched and matched groups respectively (Additional file 2: Figure S2a and Additional file 5: Table S2). This overall high genetic differentiation across the whole mitogenome was also shown by a sliding-window analysis of Fst between the two groups (Additional file 2: Figure S2d). Among the 504 differences observed, 17 were non-synonymous, giving rise to 16 amino acid changes in mitochondrial protein-coding genes (PCGs) (Additional file 5: Table S2). Using RNA-seq data, we obtained sequences of all 13 PCGs from 10 himalayanus and 4 macrurus individuals. We observed 14 fixed amino acid changes between matched and mismatched groups (Additional file 6: Table S3). Thirteen of these were shared between mismatched group and macrurus, indicating the occurrence of introgression from macrurus to himalayanus. Consistent with this, the ML tree based on sequences of concatenated 13 PCGs and PCA based on SNPs from the whole mitogenome also strongly supported the classification of five mismatched individuals with macrurus (Additional file 1: Figure S1b, Additional file 2: S2b).

Nuclear genome analysis

A further PCA based on SNPs using RNA-seq data classified all himalayanus individuals together and separated them from macrurus (Additional file 2: Figure S2c). Our sliding-window analysis of Fst across each chromosome also supported a low genetic differentiation between the two himalayanus groups with only six windows characterized by a Fst value of > 0.5 (Additional file 2: Figure S2e). Lastly, we did not detect fixed amino acid changes between the two himalayanus groups at 68 nuclear-encoded OXPHOS genes. Instead, we found eight fixed amino acid changes between himalayanus and macrurus occurring at seven genes (Additional file 7: Table S4).

Differential expression analysis

To investigate transcriptional response to mitonuclear mismatch in natural populations, we examine changes of gene expression between mitonuclear mismatched and matched groups (see Fig. 1c) across six adult tissues. Overall, we identified a small number of differentially expressed genes (DEGs) with a total of 60 DEGs across all six tissues. Among these six tissues, muscle showed the largest number with 43 DEGs, and the remaining tissues each showed fewer than 10 DEGs (Fig. 2 and Table 2). No DEGs were common to all six tissues.

Fig. 2

Volcano plots showing overall gene expression in each tissue between matched and mismatched groups with expression fold changes and adjusted P values for each gene. Differentially expressed genes that were discussed in the text were shown with their names

Table 2 Differentially expressed genes (DEGs) identified in each tissue

Functional enrichment analysis on DEGs in muscle revealed 103 significant GO terms (Additional file 8: Table S5). After reducing redundancy, we obtained 27 GO terms and most of them are related to immune response including cell killing, positive regulation of cytokine production, lymphocyte activation, leukocyte proliferation, and interferon-gamma production (Fig. 3). We also found terms associated with actin-myosin filament sliding and calcium ion transport (Fig. 3 and Additional file 8: Table S5). In contrast to muscle, we did not find any significant GO terms among DEGs identified in the other five tissues examined.

Fig. 3

Multidimensional scaling plot showing shared GO terms enriched for differentially expressed genes identified in muscle. Clustering was conducted based on semantic similarity of GO terms. The colour and size of circles correspond to the q-value calculated by Metascape and the frequency of the GO term in GO annotation database, respectively. Only the highly shared GO terms are shown with names

Discussion

Mitochondrial function relies on epistatic interactions among mitochondrial and nuclear genes. This mitonuclear interaction leads to co-adapted mitonuclear genotypes in isolated populations or lineages. In cases where allopatric populations show high levels of divergence in their mitogenomes, then mtDNA introgression on secondary contact can disrupt the co-adapted mitonuclear genotypes and lead to mitonuclear mismatches in hybrid populations [7, 8]. It is known that mitonuclear mismatch can reduce the efficiency of ATP production in the OXPHOS process and promote higher level of ROS production [38, 39], with potential negative fitness consequences [27, 59]. Yet while introgression of mtDNA has been commonly documented [60, 61], very few studies have examined the effects of this mtDNA introgression (i.e. mitonuclear mismatch) on nuclear gene expression in natural conditions.

In this study, we conducted, to our knowledge, the first assessment of the transcriptional effects of mismatched mitonuclear genotype within species across multiple adult tissues. We sampled one R. a. himalayanus population from the secondary contact region of R. a. himalayanus and R. a. macrurus, and found that half of all sampled himalayanus individuals exhibited introgression of mtDNA from macrurus (Additional file 1: Figure S1, Additional file 2: S2a and Additional file 2: S2e, Additional file 5: Table S2). This introgression of mtDNA was clear in phylogenetic analyses and in our PCA (Additional file 2: Figure S2c), with several fixed amino acid changes at mitochondrial protein coding genes (Additional file 6: Table S3) observed among the two himalayanus groups and macrurus (see also [35, 37]). Next we examine the nuclear divergence between the two himalayanus groups. The results of the PCA and sliding-window Fst analysis supported the overall similarity across the nuclear genome between these two groups (Additional file 2: Figure S2c and Additional file 2: Figure S2e) and our recent study based on a wider range of sampling also indicated that no nuclear DNA introgression had occurred between the two subspecies following mtDNA introgression [37].

Currently, more work is needed to determine the exact evolutionary forces underlying the introgression of mtDNA, alongside little or no trace of introgression of ncDNA (see a recent example in chipmunks [62]). One possible explanation is that the mitogenome confers an adaptive advantage (e.g. [63]). However, while signatures of positive selection have been detected in the mitogenome of R. affinis, demographic explanations could not be ruled out completely [37]. Alternatively, the mtDNA introgression might simply be a signal that has been left behind following hybridization, whereas ncDNA introgression was eroded via recombination following backcrossing [60]. In the future whole-genome resequencing data from multiple individuals will be needed to assess the level of nuclear introgression between the two mainland subspecies. In addition, whole-genome data can also be used to test for other mechanisms, such as demographic effects and selection against hybrids [18]. Nevertheless, the current findings provide strong support that the sampled R. a. himalayanus population includes two groups: one is a mitonuclear matched group whose mitogenome and nuclear genome are both from himalayanus (Nc-himalayanus:Mt-himalayanus), and the other is a mismatched group whose mitogenome is from macrurus but whose nuclear genome is from himalayanus (Nc-himalayanus:Mt-macrurus) (Fig. 1c).

Previously, the effects of mitonuclear mismatch on nuclear gene expression have been studied in Drosophila from RNA-seq data obtained from the whole organism [27]. In contrast, we examined transcriptional variation across multiple tissue types, and found clear tissue-specific effects of mismatched mitonuclear genotype. Specifically, mitonuclear mismatch showed a larger effect in pectoral muscle than in other tissues, with 43 DEGs identified in muscle but fewer than 10 in the other tissues, a pattern that is consistent with the fact that muscle usually requires more energy than other tissues for animals, which is likely to be especially true of volant species such as bats. An alternative explanation for why muscle shows greater number of DEGs than other tissues is that bulk RNA-seq, applied in this study, may have a larger effect on the identification of DEGs in heterogeneous tissues (e.g. brain and cochlear) than in homogeneous tissue (muscle). In the future, the application of single-cell transcriptomics approaches [64] might be informative in testing this possibility.

Although mtDNA from mismatched and matched groups shows only 14 fixed amino acid changes in protein-coding genes, our results demonstrate that this level of mtDNA difference still can likely have profound impacts on nuclear gene expression. In Drosophila, Mossman et al. [27] examined four different mtDNA haplotypes with 10 to 95 amino acid changes to assess the effects of mtDNA variation on nuclear gene expression, and identified similar numbers of DEGs in male comparisons. On the other hand, work on killifish revealed that only three nuclear genes were significantly differentially expressed between different mitochondrial genotypes in skeletal muscle [29]. However, because in this latter study the differences among mitochondrial genotypes (e.g. the number of amino acid changes or total mutations across the mitogenome) were not reported, the results cannot be compared directly to our findings, or to those of Mossman et al. [27]. Thus to draw general conclusions about the effects of different levels of divergent mitochondrial genotypes in natural populations, additional comparative transcriptomic studies of multiple tissues are needed.

We identified several significantly upregulated genes in mismatched individuals, the majority of which are associated with protection against the oxidative damage putatively caused by inefficiency of the OXPHOS pathway in the mitochondria. These genes may play essential roles for organisms in response to cellular stress in nature, and below we discuss some of these in each tissue specifically.

Most differentially expressed genes (DEGs) in muscle (37 of 43) showed upregulation in the mismatched group (Table 2). Among these, the highest expression was seen in DNASE1L3 (deoxyribonuclease 1 like 3). The loss of this gene has been linked to aberrations in the fragmentation of DNA molecules [65] and the formation of anti-DNA antibodies and autoimmunity in mice and human [66]. It has been suggested that DNASE1L3 could be activated by mitochondrial disruption [67]. We also found very high upregulation (4x) in the mismatched group of GPx3 (glutathione peroxidase 3). Its encoded protein is a member of a family of antioxidant enzymes that is critical for the antioxidant effect of the peroxisome proliferator-activated receptor γ in response to oxidative stress in skeletal muscle cells [68]. The gene HSPB6 encoded small heat shock protein HSPB6 which is most highly expressed in different types of muscle [69] and can be induced by oxidative stress [70]. Our functional enrichment analyses of DEGs in muscle revealed that a majority of DEGs are related to immune response, with other important terms associated with muscle contraction and calcium ion transport. Indeed, one term (GO:0033275) included three DEGs (ACTC1, MYL4, and TNNT2) with essential roles in muscle function, the former of which has also been shown to be overexpressed in mitochondrial myopathy due to mitochondrial respiratory chain deficiency [71]. Another important GO term (GO:0006816) includes six DEG, of which CXCL10 showed almost 32-fold higher expression in the mismatched group compared to matched group; this gene was previously shown to be important in mitochondrial dysfunction and cellular apoptosis [72].

Although DEGs from other tissues did not reveal any significant functional enrichment based on GO terms, some of these loci nevertheless have known roles in protecting against mitochondrial dysfunction. Specifically, among the two DEGs from heart that were upregulated in the mismatched group, ISG15 (interferon-stimulated gene 15) encodes a protein that is strongly associated with antiviral immune response [73] as well as regulation of mitochondrial OXPHOS and mitophagy processes during viral infection [74]. Recently, ISG15 has also been shown to play key roles in genome stability as a sensor of the DNA damage response and its expression is closely related to p53-mediated cellular processes [75, 76]. The second gene, IFI6 (interferon alpha inducible protein 6), is also one of the interferon-stimulated genes, and encodes a mitochondrial localized protein. IFI6 might be connected with cellular ATP production and OXPHOS efficiency, by directly promoting mitochondrial supercomplex assembly [77]. We also found upregulation in the brain of the mismatched group, with an almost four-fold change in CYP26A1, which encodes cytochrome P450, a protein that has been shown to suppress oxidative stress-mediated apoptosis [78]. One upregulated gene in the cochlea of the mismatched group, SLC40A1, encodes an iron transporter and can regulate cell oxidative stresses [79]. Among two genes (CA3 and NMRAL1) showing upregulation in the small intestine of the mismatched group, CA3 encodes carbonic anhydrase 3, which functions as an antioxidant and has been proposed to have protective roles against oxidative damage [80], while NMRAL1 (also called HSCARG) acts as a cellular redox sensor and can regulate DNA damage response caused by severe oxidative stress [81]. Finally, in contrast to other tissues, almost all DEGs identified in the liver are downregulated in the mismatched group, although it is unclear whether these genes function in oxidative stress.

Conclusions and future work

In our study system, we compared bats with near-identical nuclear genetic backgrounds but with contrasting histories of mtDNA introgression, and found significant and tissue-specific effects of mitonuclear mismatch on nuclear gene expression. Several genes upregulated in mismatched individuals encode proteins with known roles in responding to oxidative stress, making these potentially important candidates for future studies, including those of mitochondrial replacement therapy in human oocytes, a method used to treat mitochondrial diseases [82].

Because our current samples are all males, a remaining question that needs to be addressed is whether mitonuclear mismatches have similar effects on nuclear gene expression in females as they do in males. Previous studies aiming to answer this question have drawn different conclusions; some have indicated that mitochondrial variation have major effects on nuclear gene expression in males but little effects in females [83], a result that can be explained by the mother’s curse hypothesis, supporting a sex-specific selective sieve in mtDNA [83, 84]. However, other studies have demonstrated contradictory results, with smaller impacts of mtDNA polymorphism in males than in females [26, 85]. To address the sexual difference of mitonuclear mismatch effects on gene expression, more work needs to be performed on a wider range of animals.

When confirming the similarity of nuclear background in the two himalayanus groups, we focused on those nuclear-encoded OXPHOS genes because proteins encoded by these genes show direct interactions with mitochondrial encoded proteins. However, no fixed amino acid changes at these genes were found between the two himalayanus groups, indicating that nuclear-encoded mitochondrial genes (N-mt genes) might not have co-introgressed with mtDNA from macrurus to the mismatched group. Co-introgression of N-mt genes (mitonuclear co-introgression) has been considered as the most efficient way to overcome mitonuclear mismatches caused by mtDNA introgression [18], however, very few cases of mitonuclear co-introgression have been documented in natural populations [21, 86]. Recent studies have shown that co-introgression of N-mt genes may occur beyond those nuclear-encoded OXPHOS genes [2]. Further whole-genome resequencing data are needed to assess the extent of nuclear elements of macrurus in the genomes of mismatched individuals of himalayanus. In the future we also hope to investigate evolution of the nuclear genome under the strong selective force of foreign mtDNA by comparing genomic signature of differentiation between matched and mismatched himalayanus individuals [9, 87].

Mitonuclear mismatch is commonly detectable as mitonuclear discordance, and numerous studies have focused on understanding the processes that cause such discordance, such as incomplete lineage sorting or introgression. However, since mitonuclear discordance is likely an important intermediate state during the speciation, more work is needed to investigate the consequences of this discordance on survival and/or fitness of organism. Such studies could improve current understanding of the processes or factors that maintain mitonuclear discordance in natural populations [23]. While directly assessing the fitness consequences of mitonuclear discordance remains challenging for most wild animals (but see [10]), especially for those with long generation times, we show that studying the effects of mitonuclear genome mismatches on gene expression is more tractable (see also [26, 27, 85, 88]). Such approaches offer promising insights into the evolutionary changes and roles of natural selection in the process of adaptation to cellular stress caused by mitonuclear discordance in nature.