Skip to main content
Download PDF

Mitonuclear mismatch alters nuclear gene expression in naturally introgressed Rhinolophus bats

Frontiers in Zoology volume 18, Article number: 42 (2021) Cite this article

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
figure 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

Full size image

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 *
Full size table

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 |log2 (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
figure 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

Full size image
Table 2 Differentially expressed genes (DEGs) identified in each tissue
Full size table

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
figure 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

Full size image

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.

Availability of data and materials

All raw sequencing reads have deposited to NCBI’s Sequence Read Archive database (SRA) (BioProject accession numbers: PRJNA727985; PRJNA740060; PRJNA740148). The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. All data generated or analysed during this study are included in this published article and its supplementary information files.

Abbreviations

CPM:

Counts per million

Cytb:

Cytochrome b

DE:

Differential expression

DEGs:

Differentially expressed genes

DNASE1L3:

Deoxyribonuclease 1 like 3

GO:

Gene Ontology

GPx3:

Glutathione peroxidase 3

IFI6:

Interferon alpha inducible protein 6

ISG15:

Interferon-stimulated gene 15

Log2FC:

Log2 transformed fold change

ML:

Maximum-likelihood

MQ:

Mapping quality

mtDNA:

Mitochondrial DNA

nDNA:

Nuclear DNA

NJ:

Neighbor-Joining

N-mt genes:

Nuclear-encoded mitochondrial genes

OXPHOS:

Oxidative phosphorylation

PCA:

Principal component analysis

PCGs:

Protein-coding genes

ROS:

Reactive oxygen species

SRA:

Sequence Read Archive database

References

  1. Wolff JN, Ladoukakis ED, Enríquez JA, Dowling DK. Mitonuclear interactions: evolutionary consequences over multiple biological scales. Philos Trans R Soc Lond B Biol Sci. 2014;369(1646):20130443. https://doi.org/10.1098/rstb.2013.0443.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Hill GE. Mitonuclear ecology. Oxford University Press; 2019.

    Book  Google Scholar 

  3. Saccone C, Lanave C, De Grassi A. Metazoan OXPHOS gene families: evolutionary forces at the level of mitochondrial and nuclear genomes. BBA-Bioenergetics. 2006;1757(9–10):1171–8. https://doi.org/10.1016/j.bbabio.2006.04.021.

    Article  CAS  PubMed  Google Scholar 

  4. Hill GE. Mitonuclear compensatory coevolution. Trends Genet. 2020;36(6):403–14. https://doi.org/10.1016/j.tig.2020.03.002.

    Article  CAS  PubMed  Google Scholar 

  5. Bar-Yaacov D, Blumberg A, Mishmar D. Mitochondrial-nuclear co-evolution and its effects on OXPHOS activity and regulation. Biochim Biophys Acta. 2012;1819(9–10):1107–11. https://doi.org/10.1016/j.bbagrm.2011.10.008.

    Article  CAS  PubMed  Google Scholar 

  6. Rand DM, Haney RA, Fry AJ. Cytonuclear coevolution: the genomics of cooperation. Trends Ecol Evol. 2004;19(12):645–53. https://doi.org/10.1016/j.tree.2004.10.003.

    Article  PubMed  Google Scholar 

  7. Dowling DK, Friberg U, Lindell J. Evolutionary implications of non-neutral mitochondrial genetic variation. Trends Ecol Evol. 2008;23(10):546–54. https://doi.org/10.1016/j.tree.2008.05.011.

    Article  PubMed  Google Scholar 

  8. Barreto FS, Watson ET, Lima TG, Willett CS, Edmands S, Li W, et al. Genomic signatures of mitonuclear coevolution across populations of Tigriopus californicus. Nat Ecol Evol. 2018;2(8):1250–7. https://doi.org/10.1038/s41559-018-0588-1.

    Article  PubMed  Google Scholar 

  9. Healy TM, Burton RS. Strong selective effects of mitochondrial DNA on the nuclear genome. Proc Natl Acad Sci U S A. 2020;117(12):6616–21. https://doi.org/10.1073/pnas.1910141117.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Rank NE, Mardulyn P, Heidl SJ, Roberts KT, Zavala NA, Smiley JT, et al. Mitonuclear mismatch alters performance and reproductive success in naturally introgressed populations of a montane leaf beetle. Evolution. 2020;74(8):1724–40. https://doi.org/10.1111/evo.13962.

    Article  CAS  PubMed  Google Scholar 

  11. Burton RS, Pereira RJ, Barreto FS. Cytonuclear genomic interactions and hybrid breakdown. Annu Rev Ecol Evol Syst. 2013;44:281–302. https://doi.org/10.1146/annurev-ecolsys-110512-135758.

    Article  Google Scholar 

  12. Burton RS, Barreto FS. A disproportionate role for mtDNA in Dobzhansky–Muller incompatibilities? Mol Ecol. 2012;21(20):4942–57. https://doi.org/10.1111/mec.12006.

    Article  CAS  PubMed  Google Scholar 

  13. Ma H, Gutierrez NM, Morey R, Van Dyken C, Kang E, Hayama T, et al. Incompatibility between nuclear and mitochondrial genomes contributes to an interspecies reproductive barrier. Cell Metab. 2016;24(2):283–94. https://doi.org/10.1016/j.cmet.2016.06.012.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Dobelmann J, Alexander A, Baty JW, Gemmell NJ, Gruber MAM, Quinn O, et al. The association between mitochondrial genetic variation and reduced colony fitness in an invasive wasp. Mol Ecol. 2019;28(14):3324–38. https://doi.org/10.1111/mec.15159.

    Article  CAS  PubMed  Google Scholar 

  15. Lee HY, Chou JY, Cheong L, Chang NH, Yang SY, Leu JY. Incompatibility of nuclear and mitochondrial genomes causes hybrid sterility between two yeast species. Cell. 2008;135(6):1065–73. https://doi.org/10.1016/j.cell.2008.10.047.

    Article  CAS  PubMed  Google Scholar 

  16. Crespi B, Nosil P. Conflictual speciation: species formation via genomic conflict. Trends Ecol Evol. 2013;28(1):48–57. https://doi.org/10.1016/j.tree.2012.08.015.

    Article  PubMed  Google Scholar 

  17. Hill GE. Mitonuclear coevolution as the genesis of speciation and the mitochondrial DNA barcode gap. Ecol Evol. 2016;6(16):5831–42. https://doi.org/10.1002/ece3.2338.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Sloan DB, Havird JC, Sharbrough J. The on-again, off-again relationship between mitochondrial genomes and species boundaries. Mol Ecol. 2017;26(8):2212–36. https://doi.org/10.1111/mec.13959.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Bar-Yaacov D, Hadjivasiliou Z, Levin L, Barshad G, Zarivach R, Bouskila A, et al. Mitochondrial involvement in vertebrate speciation? The case of mitonuclear genetic divergence in chameleons. Genome Biol Evol. 2015;7(12):3322–36. https://doi.org/10.1093/gbe/evv226.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Baris TZ, Wagner DN, Dayan DI, Du X, Blier PU, Pichaud N, et al. Evolved genetic and phenotypic differences due to mitochondrial–nuclear interactions. PLoS Genet. 2017;13(3): e1006517. https://doi.org/10.1371/journal.pgen.1006517.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Morales HE, Pavlova A, Amos N, Major R, Kilian A, Greening C, et al. Concordant divergence of mitogenomes and a mitonuclear gene cluster in bird lineages inhabiting different climates. Nat Ecol Evol. 2018;2(8):1258–67. https://doi.org/10.1038/s41559-018-0606-3.

    Article  PubMed  Google Scholar 

  22. Zaidi AA, Makova KD. Investigating mitonuclear interactions in human admixed populations. Nat Ecol Evol. 2019;3(2):213–22. https://doi.org/10.1038/s41559-018-0766-1.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Streicher JW, Day JJ. The toad’s warts: discordance creates bumpy expectations of mitochondrial-nuclear evolution between species. Mol Ecol. 2020;29(18):3400–2. https://doi.org/10.1111/mec.15568.

    Article  PubMed  Google Scholar 

  24. Lee-Yaw JA, Jacobs CGC, Irwin DE. Individual performance in relation to cytonuclear discordance in a northern contact zone between long-toed salamander (Ambystoma macrodactylum) lineages. Mol Ecol. 2014;23(18):4590–602. https://doi.org/10.1111/mec.12878.

    Article  PubMed  Google Scholar 

  25. Mossman JA, Biancani LM, Zhu CT, Rand DM. Mitonuclear epistasis for development time and its modification by diet in Drosophila. Genetics. 2016;203(1):463–84. https://doi.org/10.1534/genetics.116.187286.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Mossman JA, Tross JG, Li N, Wu Z, Rand DM. Mitochondrial-nuclear interactions mediate sex-specific transcriptional profiles in Drosophila. Genetics. 2016;204(2):613–30. https://doi.org/10.1534/genetics.116.192328.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Mossman JA, Ge JY, Navarro F, Rand DM. Mitochondrial DNA fitness depends on nuclear genetic background in Drosophila. G3 Genes Genom Genet. 2019;9(4):1175–88. https://doi.org/10.1534/g3.119.400067.

    Article  CAS  Google Scholar 

  28. Flight PA, Nacci D, Champlin D, Whitehead A, Rand DM. The effects of mitochondrial genotype on hypoxic survival and gene expression in a hybrid population of the killifish, Fundulus heteroclitus. Mol Ecol. 2011;20(21):4503–20. https://doi.org/10.1111/j.1365-294X.2011.05290.x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Healy TM, Bryant HJ, Schulte PM. Mitochondrial genotype and phenotypic plasticity of gene expression in response to cold acclimation in killifish. Mol Ecol. 2017;26(3):814–30. https://doi.org/10.1111/mec.13945.

    Article  CAS  PubMed  Google Scholar 

  30. Meiklejohn CD, Holmbeck MA, Siddiq MA, Abt DN, Rand DM, Montooth KL. An incompatibility between a mitochondrial tRNA and its nuclear-encoded tRNA synthetase compromises development and fitness in Drosophila. PLoS Genet. 2013;9(1): e1003238. https://doi.org/10.1371/journal.pgen.1003238.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Barreto FS, Burton RS. Evidence for compensatory evolution of ribosomal proteins in response to rapid divergence of mitochondrial rRNA. Mol Biol Evol. 2013;30(2):310–4. https://doi.org/10.1093/molbev/mss228.

    Article  CAS  PubMed  Google Scholar 

  32. Immonen E, Rönn J, Watson C, Berger D, Arnqvist G. Complex mitonuclear interactions and metabolic costs of mating in male seed beetles. J Evol Biol. 2016;29(2):360–70. https://doi.org/10.1111/jeb.12789.

    Article  CAS  PubMed  Google Scholar 

  33. Thomas SP, Suthers RA. The physiology and energetics of bat flight. J Exp Biol. 1972;57(2):317–35. https://doi.org/10.1242/jeb.57.2.317.

    Article  Google Scholar 

  34. Mao XG, Zhu GJ, Zhang SY, Rossiter SJ. Pleistocene climatic cycling drives intra-specific diversification in the intermediate horseshoe bat (Rhinolophus affinis) in Southern China. Mol Ecol. 2010;19(13):2754–69. https://doi.org/10.1111/j.1365-294X.2010.04704.x.

    Article  CAS  PubMed  Google Scholar 

  35. Mao X, He G, Hua P, Jones G, Zhang S, Rossiter SJ. Historical introgression and the persistence of ghost alleles in the intermediate horseshoe bat (Rhinolophus affinis). Mol Ecol. 2013;22(44):1035–50. https://doi.org/10.1111/mec.12154.

    Article  CAS  PubMed  Google Scholar 

  36. Mao X, Zhu G, Zhang L, Zhang S, Rossiter SJ. Differential introgression among loci across a hybrid zone of the intermediate horseshoe bat (Rhinolophus affinis). BMC Evol Biol. 2014;14(1):1–13. https://doi.org/10.1186/1471-2148-14-154.

    Article  Google Scholar 

  37. Mao XG, Rossiter SJ. Genome-wide data reveal discordant mitonuclear introgression in the intermediate horseshoe bat (Rhinolophus affinis). Mol Phylogenet Evol. 2020;150: 106886. https://doi.org/10.1016/j.ympev.2020.106886.

    Article  PubMed  Google Scholar 

  38. Rand DM, Mossman JA, Zhu L, Biancani LM, Ge JY. Mitonuclear epistasis, genotype-by-environment interactions, and personalized genomics of complex traits in Drosophila. IUBMB Life. 2018;70(12):1275–88. https://doi.org/10.1002/iub.1954.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Cocco T, Sgobbo P, Clemente M, Lopriore B, Grattagliano I, Di Paola M, et al. Tissue-specific changes of mitochondrial functions in aged rats: effect of a long-term dietary treatment with N-acetylcysteine. Free Radical Biol Med. 2005;38(6):796–805. https://doi.org/10.1016/j.freeradbiomed.2004.11.034.

    Article  CAS  Google Scholar 

  40. Sun HJ, Chen WL, Wang JY, Zhang LB, Rossiter SJ, Mao XG. Echolocation call frequency variation in horseshoe bats: molecular basis revealed by comparative transcriptomics. P Roy Soc B-Biol Sci. 1934;2020(287):20200875. https://doi.org/10.1098/rspb.2020.0875.

    Article  Google Scholar 

  41. Mao X, Zhang J, Zhang S, Rossiter SJ. Historical male-mediated introgression in horseshoe bats revealed by multilocus DNA sequence data. Mol Ecol. 2010;19(7):1352–66. https://doi.org/10.1111/j.1365-294X.2010.04560.x.

    Article  CAS  PubMed  Google Scholar 

  42. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9. https://doi.org/10.1093/molbev/mst197.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a toolset for whole-genome association and population-based linkage analysis. Am J Hum Genet. 2007;81(3):559–75. https://doi.org/10.1086/519795.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60. https://doi.org/10.1093/bioinformatics/btp324.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9. https://doi.org/10.1093/bioinformatics/btp352.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38(6):1358–70. https://doi.org/10.1111/j.1558-5646.1984.tb05657.x.

    Article  CAS  PubMed  Google Scholar 

  49. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8. https://doi.org/10.1093/bioinformatics/btr330.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Shen YY, Liang L, Zhu ZH, Zhou WP, Irwin DM, Zhang YP. Adaptive evolution of energy metabolism genes and the origin of flight in bats. Proceed Nat Acad Sci. 2010;107(19):8666–71. https://doi.org/10.1073/pnas.0912613107.

    Article  Google Scholar 

  51. Ding Y, Chen W, Mao X. The complete mitochondrial genome of Rhinolophus affinis himalayanus. Mitochondrial DNA B. 2021;6(1):164–5. https://doi.org/10.1080/23802359.2020.1856011.

    Article  Google Scholar 

  52. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60. https://doi.org/10.1038/nmeth.3317.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30. https://doi.org/10.1093/bioinformatics/btt656.

    Article  CAS  PubMed  Google Scholar 

  54. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. https://doi.org/10.1186/s13059-014-0550-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Croux C, Filzmoserb P, Oliveirac MR. Algorithms for projection-pursuit robust principal component analysis. Chemometr Intell Lab. 2007;87(2):218–25. https://doi.org/10.1016/j.chemolab.2007.01.004.

    Article  CAS  Google Scholar 

  56. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57(1):289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.

    Article  Google Scholar 

  57. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10(1):1523. https://doi.org/10.1038/s41467-019-09234-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE. 2011;6(7): e21800. https://doi.org/10.1371/journal.pone.0021800.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Chang CC, Rodriguez J, Ross J. Mitochondrial–nuclear epistasis impacts fitness and mitochondrial physiology of interpopulation Caenorhabditis briggsae hybrids. G3 Genes Genom Genet. 2015;6(1):209–19. https://doi.org/10.1534/g3.115.022970.

    Article  CAS  Google Scholar 

  60. Toews DPL, Brelsford A. The biogeography of mitochondrial and nuclear discordance in animals. Mol Ecol. 2012;21(16):3907–30. https://doi.org/10.1111/j.1365-294X.2012.05664.x.

    Article  CAS  PubMed  Google Scholar 

  61. Hill GE. Reconciling the mitonuclear compatibility species concept with rampant mitochondrial introgression. Integr Comp Biol. 2019;59(4):912–24. https://doi.org/10.1093/icb/icz019.

    Article  CAS  PubMed  Google Scholar 

  62. Sarver BA, Herrera ND, Sneddon D, Hunter SS, Settles ML, Kronenberg Z, Demboski JR, Good JM, Sullivan J. Diversification, introgression, and rampant cytonuclear discordance in Rocky Mountains Chipmunks (Sciuridae: Tamias). Syst Biol. 2021. https://doi.org/10.1093/sysbio/syaa085.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Melo-Ferreira J, Vilela J, Fonseca MM, da Fonseca RR, Boursot P, Alves PC. The elusive nature of adaptive mitochondrial DNA evolution of an arctic lineage prone to frequent introgression. Genome Biol Evol. 2014;6(4):886–96. https://doi.org/10.1093/gbe/evu059.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Kulkarni A, Anderson AG, Merullo DP, Konopka G. Beyond bulk: a review of single cell transcriptomics methodologies and applications. Curr Opin Biotechnol. 2019;58:129–36. https://doi.org/10.1016/j.copbio.2019.03.001.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Serpas L, Chan RWY, Jiang P, Ni M, Sun K, Rashidfarrokhi A, et al. Dnase1l3 deletion causes aberrations in length and end-motif frequencies in plasma DNA. Proc Natl Acad Sci U S A. 2019;116(2):641–9. https://doi.org/10.1073/pnas.1815031116.

    Article  CAS  PubMed  Google Scholar 

  66. Sisirak V, Sally B, D’Agati V, Martinez-Ortiz W, Özçakar ZB, David J, et al. Digestion of chromatin in apoptotic cell microparticles prevents autoimmunity. Cell. 2016;166(1):88–101. https://doi.org/10.1016/j.cell.2016.05.034.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Shi G, Abbott KN, Wu W, Salter RD, Keyel PA. Dnase1L3 regulates inflammasome-dependent cytokine secretion. Front Immunol. 2017;8:522. https://doi.org/10.3389/fimmu.2017.00522.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Chung SS, Kim M, Youn BS, Lee NS, Park JW, Lee IK, et al. Glutathione peroxidase 3 mediates the antioxidant effect of peroxisome proliferator-activated receptor γ in human skeletal muscle cells. Mol Cell Biol. 2009;29(1):20–30. https://doi.org/10.1128/MCB.00544-08.

    Article  CAS  PubMed  Google Scholar 

  69. Salinthone S, Tyagi M, Gerthoffer WT. Small heat shock proteins in smooth muscle. Pharmacol Ther. 2008;119(1):44–54. https://doi.org/10.1016/j.pharmthera.2008.04.005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Wang X, Zhao T, Huang W, Wang T, Qian J, Xu M, et al. Hsp20-engineered mesenchymal stem cells are resistant to oxidative stress via enhanced activation of Akt and increased secretion of growth factors. Stem Cells. 2009;27(12):3021–31. https://doi.org/10.1002/stem.230.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Tyynismaa H, Carroll CJ, Raimundo N, Ahola-Erkkilä S, Wenz T, Ruhanen H, et al. Mitochondrial myopathy induces a starvation-like response. Hum Mol Genet. 2010;19(20):3948–58. https://doi.org/10.1093/hmg/ddq310.

    Article  CAS  PubMed  Google Scholar 

  72. Singh L, Arora SK, Bakshi DK, Majumdar S, Wig JD. Potential role of CXCL10 in the induction of cell injury and mitochondrial dysfunction. Int J Exp Pathol. 2010;91(3):210–23. https://doi.org/10.1111/j.1365-2613.2009.00697.x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Perng YC, Lenschow DJ. ISG15 in antiviral immunity and beyond. Nat Rev Microbiol. 2018;16(7):423–39. https://doi.org/10.1038/s41579-018-0020-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Baldanta S, Fernández-Escobar M, Acín-Perez R, Albert M, Camafeita E, Jorge I, et al. ISG15 governs mitochondrial function in macrophages following vaccinia virus infection. PLoS Pathog. 2017;13(10): e1006651. https://doi.org/10.1371/journal.ppat.1006651.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Tigano M, Vargas DC, Tremblay-Belzile S, Fu Y, Sfeir A. Nuclear sensing of breaks in mitochondrial DNA enhances immune surveillance. Nature. 2021;591(7850):477–81. https://doi.org/10.1038/s41586-021-03269-w.

    Article  CAS  PubMed  Google Scholar 

  76. Sandy Z, da Costa IC, Schmidt CK. More than meets the ISG15: emerging roles in the DNA damage response and beyond. Biomolecules. 2020;10(11):1557. https://doi.org/10.3390/biom10111557.

    Article  CAS  PubMed Central  Google Scholar 

  77. Liu Z, Gu S, Lu T, Wu K, Li L, Dong C, et al. IFI6 depletion inhibits esophageal squamous cell carcinoma progression through reactive oxygen species accumulation via mitochondrial dysfunction and endoplasmic reticulum stress. J Exp Clin Cancer Res. 2020;39(1):144. https://doi.org/10.1186/s13046-020-01646-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Osanai M, Sawada N, Lee GH. Oncogenic and cell survival properties of the retinoic acid metabolizing enzyme, CYP26A1. Oncogene. 2010;29(8):1135–44. https://doi.org/10.1038/onc.2009.414.

    Article  CAS  PubMed  Google Scholar 

  79. Skjørringe T, Møller LB, Moos T. Impairment of interrelated iron- and copper homeostatic mechanisms in brain contributes to the pathogenesis of neurodegenerative disorders. Front Pharmacol. 2012;25(3):169. https://doi.org/10.3389/fphar.2012.00169.

    Article  CAS  Google Scholar 

  80. Monti DM, De Simone G, Langella E, Supuran CT, Di Fiore A, Monti SM. Insights into the role of reactive sulfhydryl groups of Carbonic Anhydrase III and VII during oxidative damage. J Enzyme Inhib Med Chem. 2017;32(1):5–12. https://doi.org/10.1080/14756366.2016.1225046.

    Article  CAS  PubMed  Google Scholar 

  81. Zang W, Zheng X. Structure and functions of cellular redox sensor HSCARG/NMRAL1, a linkage among redox status, innate immunity, DNA damage response, and cancer. Free Radic Biol Med. 2020;160:768–74. https://doi.org/10.1016/j.freeradbiomed.2020.09.016.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Dobler R, Dowling DK, Morrow EH, Reinhardt K. A systematic review and meta-analysis reveals pervasive effects of germline mitochondrial replacement on components of health. Hum Reprod Update. 2018;24(5):519–34. https://doi.org/10.1093/humupd/dmy018.

    Article  CAS  PubMed  Google Scholar 

  83. Innocenti P, Morrow EH, Dowling DK. Experimental evidence supports a sex-specific selective sieve in mitochondrial genome evolution. Science. 2011;332(6031):845–8. https://doi.org/10.1126/science.1201157.

    Article  CAS  PubMed  Google Scholar 

  84. Dowling DK, Adrian RE. Challenges and prospects for testing the mother’s curse hypothesis. Integr Comp Biol. 2019;59(4):875–89. https://doi.org/10.1093/icb/icz110.

    Article  CAS  PubMed  Google Scholar 

  85. Mossman JA, Tross JG, Jourjine NA, Li N, Wu Z, Rand DM. Mitonuclear interactions mediate transcriptional responses to hypoxia in Drosophila. Mol Biol Evol. 2017;34(2):447–66. https://doi.org/10.1093/molbev/msw246.

    Article  CAS  PubMed  Google Scholar 

  86. Beck EA, Thompson AC, Sharbrough J, Brud E, Llopart A. Gene flow between Drosophila yakuba and Drosophila santomea in subunit V of cytochrome c oxidase: a potential case of cytonuclear cointrogression. Evolution. 2015;69(8):1973–86. https://doi.org/10.1111/evo.12718.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Han KL, Barreto FS. Pervasive mitonuclear coadaptation underlies fast development in interpopulation hybrids of a marine crustacean. Genome Biol Evol. 2021;13(3):evab004. https://doi.org/10.1093/gbe/evab004.

    Article  PubMed  PubMed Central  Google Scholar 

  88. Barreto FS, Pereira RJ, Burton RS. Hybrid dysfunction and physiological compensation in gene expression. Mol Biol Evol. 2015;32(3):613–22. https://doi.org/10.1093/molbev/msu321.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Sun Haijian and Wang Jiaying for assistance with the field collection and Tianli Aoxue for her help in drawing the bat picture.

Funding

This work was supported by the National Natural Science Foundation of China (No. 31570378 and No. 31630008).

Author information

Author notes

  1. Yuting Ding and Wenli Chen have contributed equally to this work.

Authors and Affiliations

  1. Institute of Estuarine and Coastal Research, East China Normal University, Shanghai, 200062, China

    Yuting Ding

  2. School of Ecological and Environmental Sciences, East China Normal University, Shanghai, 200062, China

    Wenli Chen, Qianqian Li & Xiuguang Mao

  3. School of Biological and Chemical Sciences, Queen Mary University of London, London, E1 4NS, UK

    Stephen J. Rossiter

  4. Institute of Eco-Chongming (IEC), East China Normal University, Shanghai, 200062, China

    Xiuguang Mao

Authors
  1. Yuting Ding
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Wenli Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Qianqian Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Stephen J. Rossiter
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xiuguang Mao
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

XM and SR conceived the project; YD, WC and QL analysed the data; XM wrote the draft of the manuscript with the help of YD and WC; SR provided helpful discussion throughout the project and helped edit the manuscript; XM provided the grant for this project. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Xiuguang Mao.

Ethics declarations

Ethics approval and consent to participate

Our sampling procedure was approved by the National Animal Research Authority, East China Normal University (approval ID bf20190301).

Consent for publication

No applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. Phylogenetic relationships among samples of R. a. himalayanus and R. a. macrurus based on mtDNA.

Additional file 2

. Sequence differences of mitogenome and nuclear genome between mitonuclear matched and mismatched individuals of R. a. himalayanus.

Additional file 3

. Supporting Information for “Mitonuclear mismatch alters nuclear gene expression in naturally introgressed Rhinolophus bats”.

Additional file 4: Table S1

. Detailed sequencing information of RNA-seq data and alignment rate to the reference genome for each sample.

Additional file 5: Table S2

. Detailed differences in each position of the whole mitogenome between mitonuclear matched (Him-17) and mismatched (Him-16) individuals generated by whole-genome resequencing.

Additional file 6: Table S3

. Fourteen fixed amino acid changes between matched group (Nc-himalayanus:Mt-himalayanus) and either mismatched group (Nc-himalayanus:Mt-macrurus) or macrurus at mitochondrial protein coding genes.

Additional file 7: Table S4

. Fixed amino acid changes between two himalayanus groups (Nc-himalayanus:Mt-himalayanus and Nc-himalayanus:Mt-macrurus) and macrurus at nuclear-encoded OXPHOS genes.

Additional file 8: Table S5

. Significant GO terms enriched on differentially expressed genes (DEGs) identified in pectoral muscle. GO terms named in Figure 3 are shown in bold.

Additional file 9: Table S6

. List of 77 nuclear-encoded OXPHOS genes used in Shen et al. (2010).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ding, Y., Chen, W., Li, Q. et al. Mitonuclear mismatch alters nuclear gene expression in naturally introgressed Rhinolophus bats. Front Zool 18, 42 (2021). https://doi.org/10.1186/s12983-021-00424-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12983-021-00424-x

Keywords

  • Mitonuclear discordance
  • Gene expression
  • RNA-seq
  • Introgression

Frontiers in Zoology

ISSN: 1742-9994