Ancient hybridization and mtDNA introgression behind current paternal leakage and heteroplasmy in hybrid zones

Hybridization between heterospecific individuals has been documented as playing a direct role in promoting paternal leakage and mitochondrial heteroplasmy in both natural populations and laboratory conditions, by relaxing the egg-sperm recognition mechanisms. Here, we tested the hypothesis that hybridization can lead to mtDNA heteroplasmy also indirectly via mtDNA introgression. By using a phylogenetic approach, we showed in two reproductively isolated beetle species, Ochthebius quadricollis and O. urbanelliae, that past mtDNA introgression occurred between them in sympatric populations. Then, by developing a multiplex allele-specific PCR assay, we showed the presence of heteroplasmic individuals and argue that their origin was through paternal leakage following mating between mtDNA-introgressed and pure conspecific individuals. Our results highlight that mtDNA introgression can contribute to promote paternal leakage, generating genetic novelty in a way that has been overlooked to date. Furthermore, they highlight that the frequency and distribution of mtDNA heteroplasmy can be deeply underestimated in natural populations, as i) the commonly used PCR-Sanger sequencing approach can fail to detect mitochondrial heteroplasmy, and ii) specific studies aimed at searching for it in populations where mtDNA-introgressed and pure individuals co-occur remain scarce, despite the fact that mtDNA introgression has been widely documented in several taxa and populations.

In this paper, we argue that hybridization can promote mtDNA paternal leakage and heteroplasmy also indirectly, via mtDNA introgression. Following hybridization, mtDNA introgression, defined as the permanent incorporation of genes from one species into another, has been frequently reported to such an extent that complete mtDNA replacement can also occur [20][21][22][23][24][25][26] . In many cases, however, mtDNA-introgressed and pure conspecific individuals co-exist. For example, mtDNA introgression has been observed between the sea-rock pool mosquitoes Aedes mariae and Ae. zammitii in a recently established hybrid zones along the Italian coasts of the Adriatic Sea (e.g. as a consequence of human-mediated introduction). After its introduction, Ae. mariae diffused along a transect of about 20 km, coexisting in syntopy with Ae. zammitii. Reproductive isolation between the two species is not completed and ongoing hybridization occurs which led to bidirectional mtDNA introgression (mtDNA-introgressed individuals reached 25% and 14% in Ae. mariae and Ae. zammitii, respectively) 25 .
The co-occurrence of mtDNA-introgressed and pure conspecific individuals has been observed not only when hybridization is ongoing, but also when ancient introgression occurred between taxa that are currently reproductively isolated [27][28][29] . For example, ancient introgression following hybridization has been documented in a zone of parapatry between the two European treefrogs Hyla arborea and H. intermedia 27 . Genetic analyses using both nuclear and mitochondrial diagnostic markers showed the lack of current gene exchange between the two species, while introgressed alleles were observed in both species and in all markers analysed. In particular, mtDNA of H. arborea was observed in three individuals of H. intermedia 27 .
Here, we suggest that in populations where conspecific mtDNA-introgressed and pure conspecific individuals co-exist, a mating between pure and introgressed individuals could favour paternal leakage, by relaxing the egg-sperm recognition mechanisms that prevent it, like a mating between two heterospecific individuals. If so, it can be expected that paternal leakage-driven heteroplasmy could also involve species that hybridized and introgressed in the past and now are fully reproductively isolated. We tested this expectation, using as study-system two hydrenid beetle species belonging to the Ochthebius genus.
Ochthebius quadricollis and O. urbanelliae are distributed along the Mediterranean Sea coasts and occur syntopically in the same rock pools in a sympatric area along the coasts of the Italian peninsula 30,31 . They diverged during the Plio-Pleistocene 30,32 as consequence of the climatic oscillations in the western Mediterranean region, that likely led to repeated isolation and secondary contact events between the two species 32 . Genetic analyses of sympatric populations using nuclear markers showed that hybridization and introgression occurred in the past, and then ceased due to the action of natural selection against maladaptive hybridization, which led to full reproductive isolation (speciation by reinforcement). No F1 hybrids between O. quadricollis and O. urbanelliae occur today in sympatric populations and the analysis of mating couples showed the occurrence of assortative mating. Furthermore, mating trials under laboratory conditions showed a pattern of higher premating isolation in sympatric versus allopatric populations 33 , and reproductive character displacement was found in sympatric populations of O. urbanelliae 34 .
Here, to furnish evidences to our hypothesis, we first searched for past mtDNA introgression, by using a phylogenetic approach; then, we screened O. quadricollis and O. urbanelliae individuals from sympatric and allopatric populations for the occurrence of mtDNA heteroplasmy, by developing a multiplex allele-specific PCR assay.  Table S1).

Results
Phylogenetic analyses using the Maximum Likelihood (ML) approach was performed to reconstruct the relationships among haplotypes (Fig. 3). All haplotypes of the clade I were found in O. quadricollis individuals with the exception of the haplotypes h17 that was found in both O. quadricollis and O. urbanelliae individuals, and with the exception of the haplotype h35 that was found in one O. urbanelliae individual (Fig. 3, Table 1). All haplotypes of the clade II were found in O. urbanelliae individuals with the exception of the h34 haplotype that was found in both O. quadricollis and O. urbanelliae individuals (Fig. 3, Table 1).  Fig. S2).

Heteroplasmy.
The MAS-PCR assay was then used to analyze all individuals that were sequenced for the COI mtDNA gene fragment to check for mitochondrial heteroplasmy. In O. urbanelliae, two heteroplasmic individuals were found (i.e. one male and one female) ( Table 1). They showed two electrophoretic bands following MAS-PCR ( Supplementary Fig. S3), which sequences were identical to a fragment of the O. quadricollis and O. urbanelliae haplotypes. In O. quadricollis, no heteroplasmic individuals were found.

Discussion
Paternal leakage and mtDNA heteroplasmy have been documented in several animal species, including insects, fishes, reptiles, birds and mammals 8,10 . Being more difficult to detect when mtDNA variation within or among populations is low, these phenomena have been often observed in secondary contact zones or laboratory crosses, as a consequence of hybridization between heterospecific individuals or between homospecific individuals belonging to highly divergent genetic lineages 8,10,35 . Here, we hypothesized that hybridization can lead to mtDNA paternal leakage and heteroplasmy also via mtDNA introgression. By using genetic data, we showed that: i) mtDNA introgression occurred between the two species and that mtDNA-introgressed and pure individuals co-occur in natural populations; ii) heteroplasmic individuals are present in O. urbanelliae; iii) mtDNA heteroplasmy originated through paternal leakage.
First, the phylogenetic analyses showed that O. quadricollis and O. urbanelliae share mtDNA haplotypes, as well as having some haplotypes more closely related to those of the other species than to those of its own haplogroup (Fig. 3, Table 1). This pattern is consistent with a process of mtDNA introgression and it is concordant with the signatures in the nuclear genomes of introgressive hybridization between O. quadricollis and O. urbanelliae, before the completion of the speciation process by reinforcement, shown by nuclear markers 31,33 . Furthermore, alternative hypotheses, such as incomplete lineage sorting can be excluded 20 . Contrary to what would be expected under this hypothesis, the introgressed haplotypes between O. quadricollis and O. urbanelliae are indeed not randomly distributed across the species' range, but rather they are geographically localized within the sympatric area between the two species, where both mtDNA-introgressed and pure individuals co-occur (Fig. 1 Supplementary Fig. S3). Misleading results due to unspecific PCR amplification of the mtDNA of the two species can be confidently ruled out, as the MAS-PCR assay proved to be highly specific. Furthermore, the sequences obtained from the amplicons were identical to the COI sequences of O. quadricollis and O. urbanelliae found in the sympatric area, which also exclude the possibility that non-functional nuclear copies of mitochondrial genes (NUMTs) have been analysed.
Third, mtDNA heteroplasmy originated by paternal leakage. Empirical evidence showed that different routes can lead to mtDNA heteroplasmic individuals 8,10 . When mtDNA heteroplasmy is due to the leakage of paternal mitochondrial genome, it is expected that the mtDNA haplotypes within the heteroplasmic individuals derive from both maternal and paternal mtDNA lineages 8,10,14 . The heteroplasmic individuals found in O. urbanelliae are concordant with this expectation, as they showed one haplotype identical to the haplotypes of O. quadricollis, and the other identical to those of O. urbanelliae, which supports the hypothesis that they originated Taking into account the above results and the reproductive isolation among the two species, we argue that mtDNA introgression promoted paternal leakage in the Ochthebius beetles. Events of paternal leakage and heteroplasmy have often been ascribed to a failure of the recognition mechanisms between nuclear and mitochondrial genomes located in the eggs and sperms, and often they have been shown to involve hybridization between divergent lineages 8,36 . However, hybridization between O. quadricollis and O. urbanelliae cannot be invoked to explain paternal leakage and mtDNA heteroplasmy in these species, because reproductive isolation between them is completed 31,33 . Likewise, the hypothesis that paternal leakage originated from past hybridization events between O. quadricollis and O. urbanelliae and then was followed by maternal transmission of heteroplasmy can be excluded as well, because the heteroplasmic variants would have disappeared over time since the completion of reproductive isolation between the two species 30,[37][38][39] . Paternal leakage and mtDNA heteroplasmy in the O. quadricollis and O. urbanelliae beetles would therefore be due to mates between mtDNA-introgressed and pure conspecific individuals.
The role of the mtDNA introgression as a source of evolutionary novelty has become increasingly recognized in recent years [40][41][42][43] . Our results, by showing that mtDNA introgression can promote paternal leakage, highlight that introgression would be able to generate genetic novelty in a way that has been neglected to date. Furthermore, they support that the frequency and distribution of paternal leakage and mtDNA heteroplasmy can be deeply underestimated in natural populations, as the commonly used PCR-Sanger sequencing approach can fail to detect mitochondrial heteroplasmy, and specific studies aimed at searching for heteroplasmy in populations where mtDNA-introgressed and pure individuals co-occur remain scarce. www.nature.com/scientificreports www.nature.com/scientificreports/ In animals, the mtDNA genome has recently been shown also to recombine, both at intra-and intermolecular level 8,14,44 , and recombination has been suggested as an important evolutionary mechanism avoiding deleterious mutation meltdown in mtDNA (i.e. the Muller's ratchet), as well as for the origin of new genetic combinations. Paternal leakage and mtDNA heteroplasmy, allowing two different genomes to meet, represent an intermediate but critical step for recombination. An understanding of how frequently they occur in nature, as well as the processes underlying them, represent therefore outstanding questions to be addressed.

Materials and Methods
Samples and species recognition. Ochthebius quadricollis and O. urbanelliae adult individuals were collected from 16 localities along the Italian coasts (Table 1, Fig. 1). All individuals were identified using the morphological and biochemical keys of Audisio et al. 45 and Urbanelli et al. 30 . Standard horizontal starch gel electrophoresis was performed following the protocols described in Urbanelli 31 . mtDNA sequencing and phylogenetic inferences. DNA from single adults was extracted from the tissue homogenate used for allozymic analyses by standard CTAB (cetyltrimethyl ammonium bromide) protocol 46 . Partial sequences of the Cytochrome Oxidase I (COI) mitochondrial gene were obtained through PCR-amplification. The universal primers pair C1-J-2183 and TL2-N-3014 was used to amplify and sequencing a COI fragment of ~900 base pairs 47,48 . The following specific primers were then designed and used for further analyses: OchtCOI-f 5′-accaggatttggaataattt-3′ and OchtCOI-r 5′-tccaatagaagaaataatatttc-3′. The PCR was carried out in a 25 µl volume containing 5 ng of DNA, 10 mMTris-HCl, pH 8.3, 2.0 mM MgCl 2 , 0.2 μM of the forward and reverse primers, 0.4 mM dNTPs, 0.3 units of NZYTaq polymerase (NZYtech, Lisbon, Portugal) and water. Negative controls including all reagents but DNA were also included in the reaction. PCR cycling www.nature.com/scientificreports www.nature.com/scientificreports/ procedure was: 95 °C for 5 min followed by 34 cycles of 93 °C for 1 min, 55 °C for 1 min, 72 °C for 1 min 30 s, and a single final step at 72 °C for 10 min. PCR sequences were obtained using an ABI PRISM 3700 DNA sequencer by GATC Inc. (www.GATC.com). All individuals were double strand sequenced. The individuals found introgressed were amplified multiple times and two independent PCR products were sequenced.
Sequences were edited and aligned using the software Chromas 2.6.5 (Technelysium, Helensvale, Australia) and ClustalX 2.1 49 , respectively, and compared with the O. quadricollis and O. urbanelliae sequences available in GenBank using BLASTN algorithm. Nucleotide and amino-acidic polymorphisms of the COI gene fragments were estimated using the software DNAsp 6 50 . The average uncorrected p-distance between groups of haplotypes and populations was computed using Mega 7.0 51 .
Phylogenetic relationships among haplotypes of O. quadricollis and O. urbanelliae were inferred using the Maximum Likelihood (ML) method, as implemented in PAUP 4.0 52 . ML analysis was performed using heuristic searches with 100 rounds of random sequence addition and tree bisection-reconnection (TBR) branch swapping algorithm. The GTR + G + I substitution model was used (gamma shape parameter G = 1.480; proportion of invariable sites I = 0.669), following the substitution model inferred by jModelTest 2.0 53 using the Akaike Information Criterion. The robustness of the inferred ML tree topology was assessed by the non-parametric bootstrap method with 1000 replicates. Ochthebius heeri heeri (KT804242.1) and Ochthebius metallescens (HF931191.1) were used as outgroups in the phylogenetic analyses 54 .  (Table 1); the haplotype h35 (in bold) was found in one O. urbanelliae individual (Table 1), but is more closely related to the O. quadricollis haplotypes.