Genomic consequences of human‐mediated translocations in margin populations of an endangered amphibian

Abstract Due to their isolated and often fragmented nature, range margin populations are especially vulnerable to rapid environmental change. To maintain genetic diversity and adaptive potential, gene flow from disjunct populations might therefore be crucial to their survival. Translocations are often proposed as a mitigation strategy to increase genetic diversity in threatened populations. However, this also includes the risk of losing locally adapted alleles through genetic swamping. Human‐mediated translocations of southern lineage specimens into northern German populations of the endangered European fire‐bellied toad (Bombina bombina) provide an unexpected experimental set‐up to test the genetic consequences of an intraspecific introgression from central population individuals into populations at the species range margin. Here, we utilize complete mitochondrial genomes and transcriptome nuclear data to reveal the full genetic extent of this translocation and the consequences it may have for these populations. We uncover signs of introgression in four out of the five northern populations investigated, including a number of introgressed alleles ubiquitous in all recipient populations, suggesting a possible adaptive advantage. Introgressed alleles dominate at the MTCH2 locus, associated with obesity/fat tissue in humans, and the DSP locus, essential for the proper development of epidermal skin in amphibians. Furthermore, we found loci where local alleles were retained in the introgressed populations, suggesting their relevance for local adaptation. Finally, comparisons of genetic diversity between introgressed and nonintrogressed northern German populations revealed an increase in genetic diversity in all German individuals belonging to introgressed populations, supporting the idea of a beneficial transfer of genetic variation from Austria into North Germany.


| INTRODUC TI ON
Genetic admixture between populations is a widespread phenomenon that instantly increases the diversity of a population's gene pool and can enable a rapid adaptation to change in environmental conditions (Arnold & Kunte, 2017;Janes & Hamilton, 2017). However, in species with low dispersal abilities, natural gene flow may be limited among fragmented and isolated populations, especially those situated at the species' range margin. When facing rapid environmental changes, the adaptive potential at the margins may hence be compromised, as these genetically depauperate populations may suffer from genetic erosion and inbreeding (Arrendal et al., 2004;Madsen et al., 1999;Mussmann et al., 2017;Westemeier et al., 1998). Here, "genetic rescue" can be implemented as a mitigation strategy through translocation of individuals between populations (Johnson et al., 2010;Vilà et al., 2003). The transmission of beneficial gene variants from a closely related species or distant population is called "adaptive introgression" (Arnold & Kunte, 2017;Hedrick, 2013). Adaptive introgression could be particularly important in response to climate change as relevant alleles required at a species' poleward margin may already be present in more central populations. Introducing gene variants adapted to warmer conditions could potentially enhance the ability of the recipient poleward population to adapt more rapidly (Bridle & Vines, 2007).
One species with range margin populations that are particularly vulnerable is the fire-bellied toad (Bombina bombina), a small amphibian species found in the Central and Eastern European lowlands. It mostly inhabits permanent, fish-free freshwater bodies, like shallow stagnant lakes, ponds and swamps (Drews et al., 2011). After the Pleistocene, this species likely colonized Europe in two distinct expansions out of refugia located north-west of the Black Sea (Fijarczyk et al., 2011;Hofman et al., 2007). The Carpathian Mountains and the Central European middle mountains at the border between Germany and the Czech Republic formed natural barriers that led to a north-south division of the species range. As a consequence, two genetically distinct evolutionary lineages of B. bombina emerged (Hofman et al., 2007).
Within the respective lineages, even geographically close populations of B. bombina are often genetically distinct, likely due to high site fidelity and limited dispersal capability Engel, 1996;Schröder et al., 2012). This is particularly true for locally restricted and highly fragmented populations at the northwestern margin of the species range, which runs from the south of the Swedish province Scania to the Baltic Sea through Denmark and Germany (Schleswig-Holstein, Niedersachsen, Sachsen-Anhalt) . In comparison with their counterparts in the central distribution range, margin populations of B. bombina are severely threatened by habitat alterations, pollution and intensive agriculture (Drews et al., 2011). Continuous human-mediated habitat fragmentation has led to a reduction in gene flow and locally isolated populations are genetically depauperate .
B. bombina population sizes in Germany have been declining since the middle of the last century, which led to the species conservation status being listed as "critically endangered" within Germany (Kühnel et al., 2009).
Genetic analyses from toads collected in northern Germany in 2006 uncovered an introgression of southern lineage (putatively Austrian) genotypes into northern lineage populations. This was likely the result of occasional (and illegal) translocations of toads of southern origin into northern populations .
This study was, however, restricted to a short fragment of the mitochondrial Control Region and putatively neutral nuclear microsatellites, precluding inferences of the adaptive potential of these introgression events.
Here, we test the genetic impact of this translocation on the recipient populations on a genomic scale. Specifically, we searched entire mitochondrial genomes and nuclear transcriptomes for signs of introgression from southern lineage B. bombina from Austria into conspecific northern German B. bombina populations. We inferred levels of introgression in the nuclear genome and evaluated the resultant adaptive consequences by identifying introgressed alleles spreading in the northern German gene pool, as well as local northern German alleles retained, putatively due to local adaptation. We further investigated transcriptome-wide genetic diversity of both introgressed and autochthonous populations.

| Samples
Thirty wild tadpoles were collected, five from Austria (Lobau, German tadpoles were frozen in liquid nitrogen and stored at −80°C. The five Austrian tadpoles were transferred into a RNAlater storage solution (Sigma-Aldrich Chemie GmbH).

| RNA extraction, mRNA enrichment and cDNA library preparation
We extracted RNA from the tail muscle of the fins, using a modified Trizol isolation protocol in combination with a column-based commercial kit (QIAGEN RNeasy Mini Kit) based on the manufacturer's protocol. Only skin and muscle tissue of the tadpole's fins were used in order to prevent contamination from gut bacteria during sequencing as much as possible. To further increase RNA yield, we transferred the frozen tadpole tails into a 1 ml Trizol solution and performed a homogenization with a Tissue Lyzer (6 min, 50 HZ).
We implemented the column-based clean-up procedure from the QIAGEN Rneasy Mini Kit as described in the manufacturer's manual.
We assessed the quality and quantity of the extracted total RNA on an Agilent Bioanalyzer 2100. We only selected samples with a high RIN number (>8) as well as a high RNA concentration (c > 50 ng/µl) for library preparation. We then performed an mRNA enrichment on these samples with poly(A) beads using the Bio Scientific NEXTflex

| Quality filtering and nuclear mapping
We trimmed Illumina adapter sequences and poly A tails, only keeping reads with a quality score of at least 25 and a minimum length of 30 bp for all 30 samples using the software Skewer v0.2.2 (Jiang et al., 2014). We mapped the trimmed reads twice using Burrows-Wheeler algorithm (BWA) v0.7.15 (Li & Durbin, 2009) specifying the mem algorithm and default parameters to a previously published reference transcriptome of B. bombina from NCBI (GenBank accession: HADQ00000000.1) for population genomics analyses, and a transcriptome from B. orientalis (GenBank accession: HADT00000000.1) for use in downstream admixture analyses to avoid ascertainment bias that can be caused by mapping to an ingroup reference , as well as a transcriptome from Bombina variegata (GenBank accession: HADR00000000.1).

| Mitogenomes
We additionally mapped the processed reads using BWA twice independently to two published reference B. bombina mitogenomes, one from Austria (GenBank accession: JX893173.1) (Pabijan et al., 2013) and the other from Germany (GenBank accession: MH893761.1) (De Cahsan et al., 2019). Due to insufficient and low read coverage, we removed one sample from the German Högsdorf population (Ge-Ho-KQ01). Mitochondrial consensus sequences were built using a maximum effective base depth approach in ANGSD (-doFasta 3) (Korneliussen et al., 2014) and specifying the following parameters -minMapQ 25 -minQ 25 -uniqueOnly 1. We then compared the total number of mapped reads per individual for each reference (Table   S1). In order to rule out ascertainment biases, which could be caused by mapping all individuals to a single reference from an ingroup population, we additionally constructed two independent maximumlikelihood phylogenetic trees, one built using all individuals mapped to the Austrian reference and one using all individuals mapped to the German reference ( Figure S2). This was done by computing an alignment of our 29 B. bombina samples and two outgroup mitogenome sequences (B. orientalis (GenBank accession: AY957562.1), B. variegata (GenBank accession: AY971143.1)) using MAFFT v7.392 (Katoh & Standley, 2013) and running RAxML v8.2.11 (Stamatakis, 2006(Stamatakis, , 2014, with 1000 bootstrap iterations, specifying a GTR+G substitution model. The resultant trees were then visually inspected for incongruencies. As the resultant trees were the same regardless of reference, we used the reference that had most reads mapping to it to construct a consensus sequence for each individual. These consensus sequences were aligned using MAFFT. We then produced an unrooted maximum-likelihood phylogenetic tree from these aligned consensus sequences by running RAxML with 1000 bootstrap iterations, specifying a GTR+G substitution model.

| Population transcriptomic analyses
We performed two independent principal component analyses with the software ANGSD v0.921 and PCangsd v0.95 (Meisner & Albrechtsen, 2018), one for all B. bombina individuals and one including only individuals from Germany. For both analyses, we used a minimum base quality (-minQ) and a minimum mapping quality (-minMapQ) parameter of 25; we required 20 individuals as the minimum number for a site to be considered for analysis , only included reads mapping uniquely to one location (-uniqueonly), specified a minimum minor allele frequency of 0.05 (-minFreq 0.05), called genotype likelihoods using the GATK method (-GL 2), only included contigs over 1000 bp in size, required each site to have at least five reads mapping (-mininddepth 5) and only considered SNPs with a p-value lower than 1e-6 (-SNP_pval). The output genotype likelihoods were then converted from a beagle file into a covariance matrix using PCangsd. The same genotype likelihoods were also used to perform an admixture proportions analysis in PCangsd (-admix).
The most probable K value was chosen by the software based on the eigenvalues of the PCA.

| Transcriptome-wide D-statistics/ABBA-BABA analysis
We ran D-statistics in ANGSD to look for the relative degree of gene flow in our samples. D-statistics works based on a predefined species tree which uses three ingroup and one outgroup taxon.
This topology can be written as (((H1,H2),H3),O) where H1 and H2 are more closely related to one another than either are to H3. Dstatistics then scans across the genome to find regions that break the known species tree, either due to incomplete lineage sorting or admixture between H3 and either H1 or H2. The occurrence of (((H1,H3),H2),O) being equal to (((H3,H2),H1),O), is interpreted as incomplete lineage sorting (or equal amounts of gene flow between H1+H3 and H2+H3). Sites at which H2 and H3 share a derived allele, while H1 and the outgroup share the ancestral allele, are called ABBA sites, while BABA represents sites at which H1 and H3 share the derived state and H2 and the outgroup share the ancestral allele.
Any deviations from the expected 1:1 ratio are considered as differential levels of gene flow between the ingroup species analysed and are represented by the D value (null expectation of the test: D = 0).
This value quantifies the deviation from the expected ratio, which is calculated by the differences in the sum of ABBA and BABA site patterns across the genome, divided by their sum: We used a random base sampling approach (-doabbababa 1), considered only transcripts greater than 1000 bp in length and used the following filters: -minMapQ 25 -minQ 25 -uniqueOnly 1. As the sister species B. variegata is known to occasionally hybridize with B. bombina in our sampling area in north-eastern Austria (Gollmann, 1984), we specified the related East Asian B. orientalis as outgroup.
Because B. orientalis is more distantly related, it fits less well to the expectations of the infinite site model assumed by ABBA/BABA.
The D-statistics has, however, been found to be robust to small deviations from the infinite site model regarding the outgroup (Zheng & Janke, 2018). To investigate the significance of our result, we performed a weighted block jackknife test specifying a block size of a single contig. As is most commonly practised with determining the significance of D-statistics, D values more than three standard errors different from zero (−3 < Z > 3) were considered as statistically significant, which corresponds to a p-value <0.0013 (Green et al., 2010;Patterson et al., 2012;Zheng & Janke, 2018).

| Transcriptome-wide Treemix analysis
To further investigate the presence of admixture between Austrian and German lineages, we ran Treemix v1.1.3 (Pickrell & Pritchard, 2012) on the transcriptome data. To produce the input file for treemix, we first called pseudo haploid sequences using ANGSD, using the same parameters as for the PCA above with the addition of -doHaploCall 2 to call the consensus haploid base. This was done in a single run specifying all individuals from our study and the two outgroup individuals from NCBI (B. variegata, and B. orientalis). The output haploid file was converted to PLINK format using haplotoplink from the ANGSD toolsuite. We used plink v1.90b3.42 (Purcell et al., 2007) to convert the ANGSD output plink file into a Freq file, which was then converted into a treemix input file using plink2treemix.py (https://bitbu cket.org/nygcr esear ch/treemix). We first ran treemix specifying no migration edges and B. orientalis as root to produce a consensus tree. We then used this tree as a prior when specifying multiple migration edges (-m 1-4) while running treemix. The tree and residues were produced in R using the plotting_funcs.R script available with treemix (Figure 2a,b, and S5a-h).

| Genetic diversity
We estimated two independent measures of genetic diversity in our individuals, observed heterozygosity and nucleotide diversity. Heterozygosity gives an indication of the individual levels of diversity while the nucleotide diversity gives an indication of the population-level diversity. We calculated heterozygosity using site allele frequencies determined using ANGSD. This was done independently for each individual using the same parameters as for the PCA with the addition of -doSaf 1. The resultant file was then run through realSFS from the ANGSD toolsuite with the parameter -nSites 50,000 to output the number of homozygous and heterozygous sites in windows of 50 kb of covered bases, traversing contig boundaries. To obtain population-level estimates, we concatenated the results for all windows from all individuals within a single population.
For the nucleotide diversity estimates, we used the same pseudo haploid call file from the previous section. We then converted the output haploid call file into a geno file to be run through a suite of python scripts available online (https://github.com/simon hmart in/ genom ics_general). While running the popgenWindows.py, we specified the window size as 20 kb so that the window size encompasses complete contigs/genes, and a minimum number of sites per window as 100. The significance of differences between populations was calculated using a Wilcoxon-Mann-Whitney test while implementing a Bonferroni correction (Table S4).

| Alleles highly differentiated between German and Austrian populations
The popgenWindows.py python script simultaneously calculated population differentiation statistics (F st ) resulting in independent F st values for each contig/transcript. From this output, we found contigs that were highly differentiated between the German and Austrian populations, in an attempt to uncover German/Austrian-specific genes. We took the contigs showing the highest F st values (top 5%) for each pairwise comparison of a German population against the Austrian population. This resulted in five pairwise comparisons (i.e. Austria-Fehmarn, Austria-Eutin, Austria-Dannau, Austria-Testorf and Austria-Högsdorf). We extracted the contigs that were found to be in the top 5% of all pairwise comparisons and cross-referenced the sequences using BLAST to identify names and functions of these genes (Table S2 and S3).

| Significantly introgressed alleles in the transcriptome
To find genes that exhibit a significant elevation of introgressed alleles in German populations, we ran the ABBABABAwindows.py (Martin et al., 2015) on the geno file from above, specifying B. orientalis as outgroup, all individuals from Austria as the H3 position, all individuals from had been reported, physically separated from the other populations by the surrounding Baltic Sea. We specified the window size as 20 kb so that the window size encompasses complete contigs/genes, and a minimum number of sites per window as 100. This method computes the F d statistic which is more robust to biases caused by genetic differentiation and nucleotide diversity when estimating D with small window sizes (Martin et al., 2015). We extracted the genes showing the strongest signs of admixture relative to Fehmarn (top 5%) and compared them between populations. From this, we estimated whether highly significant signs of admixture were found in the same contigs in all populations or whether they differed between populations, with prevalence of particular introgressed contigs in all introgressed populations putatively indicating a selective advantage. We used the online software tool panther (Protein ANalysis THrough Evolutionary Relationships) Classification System v14.1 ( Mi et al., 2013;Thomas et al., 2003) to assign genes with high levels of introgression to Gene Ontology classes.

| Nuclear population structure
To examine whether the pattern found in the mitochondrial genomes is also reflected in the coding regions of the nuclear genome, the same cDNA libraries that we sequenced and used to reconstruct mitochondrial trees were also used to build transcriptomes for each individual. We ran a transcriptome-wide principal component analy-  (Figures 1 and S3), indicating that these specimens do not originate from a common gene pool in Hardy-Weinberg equilibrium, but rather each carries unique combinations of introgressed versus local alleles throughout their genomes.

| Admixture analyses
In combination, our findings based on the mitochondrial genomes and the PCA on nuclear data suggest admixture between Austrian and some German populations. To identify the degree of admixture between southern (Austria) and northern (Germany) European fire-bellied toads, we performed an admixture structure analysis using PCangsd (Meisner & Albrechtsen, 2018). Our treemix analyses also confirm the result that the Högsdorf population is admixed with Austrian toads. When specifying no migration edges, we obtained the expected tree topology with all German populations being sister to the Austrian population. The residue matrix for this tree had a high SE between Högsdorf and Austria indicating that a migration edge between the two may be missing. We therefore ran treemix with 1, 2, 3 and 4 migration edges. When specifying a single migration edge, as expected due to the residue matrix of m = 0, we uncover migration between Austria and Högsdorf (Figure 2a,b).
However, when specifying more than one migration edge, the tree topology changes and Austria and Högsdorf become a sister clade to the remaining German populations. This is likely due to the high admixture proportion of Austrian alleles in Högsdorf.

| Genetic diversity
We independently calculated observed heterozygosity for all B. bombina individuals in this study (Figure 3a). We found the highest average level of heterozygosity in individuals from the German population of Högsdorf and in the Austrian toads. The nucleotide diversity (pi) follows a similar pattern (Figure 3b). When considering both heterozygosity and pi, we found the difference in diversity levels between Högsdorf and Austria to be nonsignificant, while these two populations had significant differences to the remaining four German populations (Table S4).

| Highly differentiated genes between German and Austrian populations
Using the population differentiation statistics (F st ), we were able to identify loci that were highly differentiated (top 5%) between the German and Austrian populations. We uncovered 23 contigs found to be highly differentiated in all pairwise comparisons between the German and Austrian populations. We cross-referenced these contigs via BLAST and identified gene names for 17 of the 23 genes of interest (Table S2). All 17 identified loci belong to genes that appear to play an important role for the general functioning of cells. However, one gene (FAM3A) is additionally associated with antifungal properties.

| Genes showing high levels of introgression
To uncover which introgressed alleles may be adaptive, we utilized F d statistics (Martin et al., 2015). We uncovered a number of genes showing the strongest signs of admixture (top 5% F d ) relative to Fehmarn (considered autochthonous, see above) in the remaining German populations (Figure 4). We identified 21 loci, showing consistent strong signals of admixture in all four German populations, excluding Fehmarn. We successfully found matches for eight out of the 21 loci using BLAST and were able to classify five of them (Table S3 and Figure S6). Interestingly, two of these genes are known to be involved in immune and stress response-related processes, as well as metabolic pathways. The

| D ISCUSS I ON
The question of whether the benefits of translocations as a method for genetic rescue outweigh the downfalls is a highly relevant topic as many species now require at least some means of human-mediated intervention for their survival (Barton, 2001;Mansfield & Land, 2002;Seddon, 2010). We sought to address this question through the production of complete mitochondrial genomes and transcriptomes from an endangered amphibian species, the European fire-bellied toad. For this, we sampled five endangered north- We confirmed the presence of southern "alien" (i.e. Austrian) DNA in the northern populations as was previously found using toads collected in 2006 . However, in contrast to the previous study where the introgression was found limited to Another scenario is migration of introgressed B. bombina individuals between sites over the last 10 years, spreading the introgressed alleles to neighbouring populations. The latter appears probable as adult toads regularly migrate up to 500 m and occasionally more than 1 km over a period of multiple days (Engel, 1996;Günther & Schneeweiss, 1996). Closely related B. variegata have been shown to migrate up to 90 m per night (Abbühl & Durrer, 1996).
Next, to uncover whether this introgression may have had adaptive consequences, we investigated which specific Austrian alleles were introgressed into the northern German populations and how widespread they became among the German populations. We found a number of gene variants present in all of the four admixed German populations indicating that they may have some selective advantage to the recipient toads leading to their retention ( Figure 4 and Table   S3). Alternatively, these gene variants could have become locally abundant due to random effects (i.e. genetic drift) and subsequently spread by dispersal. However, in favour of a potential selective advantage of the introgression event(s), two of these genes (MTCH2, and DSP) display functions that could putatively be linked to differential fitnesses. The MTCH2 gene is highly expressed in white adipose tissue (Kulyté et al., 2011;Speliotes et al., 2010) and has been linked to obesity in humans through genomewide association and expression studies (Kulyté et al., 2011;Speliotes et al., 2010).
The Austrian MTCH2 gene variant may represent an adaptive advantage for northern German B. bombina toads as it could influence the storage efficiency of fat as a form of energy. In amphibians, a greater body mass often correlates with higher reproductive success (Feder & Burggren, 1992;Kuramoto, 1978;MacCracken & Stebbings, 2012).
DSP encodes a unique desmosomal protein vital for epidermal development (Garrod & Chidgey, 2008) and has been shown to play a critical role in Xenopus laevis embryos during the frog's epidermal morphogenesis (Bharathan & Dickinson, 2019). DSP appears to be required to move basally located cells into the outer epidermal layer, where they differentiate into specific epidermal cell types and are fundamental for the skin's integrity (Bharathan & Dickinson, 2019).
An intact and resilient skin is vital for amphibians, especially at the tadpole stage as it protects them from pathogens through secretion as well as from mechanical stress. A reduction or lack of DSP in X. laevis embryos led to skin fragility and reduced resistance to mechanical stress (Bharathan & Dickinson, 2019). Although we did not detect amino acid changes within these two genes consistent with gene flow between northern and southern B. bombina populations, upstream cis-regulatory elements were likely also introgressed due to linkage disequilibrium which may influence function. However, we are unable to look for lineage-specific cis-regulatory elements using the dataset at hand.
Despite the presence of these putatively beneficial introgressed alleles, the question remains in how far genetic swamping may have compromised local adaptation in these introgressed populations, putatively reducing the fitness of admixed individuals. Therefore, we also searched for genes where local northern German alleles may have been retained despite admixture, potentially due to a role in adaptation to the local environment. We found that most such genes play an essential role for the general functioning of cells and therefore might be crucial for the German toads. However, one gene (FAM3A) stood out as it encodes for a peptide known as Drosomycin-Like Defensin (DLD). Studies on humans and Drosophila melanogaster discovered that DLD may act as a defence against invading fungal microorganisms (Simon et al., 2008). It seems reasonable to assume that this antifungal activity may be a sign of German toads having evolved a defence mechanism against the local fungal pathogens, as this allele remains present in the northern German gene pool despite ongoing admixture. This is in line with the finding of a previous study that local alleles are retained in genes of the adaptive immune system (MHC), despite introgression .
The presence of introgressed potentially beneficial alleles at some genes, but local alleles at other genes in the recipient northern German populations could constitute the result of natural selection, leading to mosaic genotypes with favourable alleles at multiple loci, be they introgressed or autochthonous. We further found that the introduction of allochthonous individuals led to an increase in genetic diversity of the recipient populations ( Figure 3). This increase in genetic diversity may provide a new lifeline for the endangered range margin populations of this species as they are often highly fragmented, have low levels of genetic diversity, and lack gene flow from surrounding populations, Although the translocation of individuals between populations can be risky to recipient populations, by introducing foreign parasites or pathogens (Cunningham, 1996;Kock et al., 2010) or disrupting locally adapted gene assemblages via an influx of foreign alleles (Edmands, 2007), it can also be beneficial. Adaptive introgression can instantly increase the genetic diversity of depauperate populations by providing new adaptive gene variants and reducing the level of inbreeding and therefore the chance of rare deleterious alleles being present in a homozygous state (Charlesworth & Willis, 2009;Edmands, 2007) which may be the case in fragmented northern fire-bellied toad populations. We show that the translocation of individuals from Austria had a significant genetic impact on northern German populations, with a tendency of spreading further. We find evidence that the introduction of alleles from Austria into Germany may have been beneficial to the recipient populations through the increase in genetic diversity and introduction of new alleles that were subsequently retained at high frequencies. Moreover, positive population trends have been reported for these introgressed populations (Drews et al., 2011), further corroborating our findings and our hypothesis that the introgression of southern alleles had a positive effect on the endangered northern populations.

| CON CLUS IONS
Through the production of mitochondrial genomes and transcriptomic data from six populations of B. bombina, we were able to assess the genetic consequences of a human-mediated translocation of individuals from the more genetically robust central Austrian population into endangered range margin populations in northern Germany.
Although conservationists may be reluctant to translocating individuals from different evolutionary lineages as means of genetic rescue (Kock et al., 2010;Schröder et al., 2012), we show that this introgression may have been beneficial to the recipient populations by introducing new alleles at some loci, while maintaining locally adaptive alleles at others, overall increasing the genetic diversity of the population. Unless local populations are completely swamped by allochthonous individuals, crossbreeding between allochthonous and autochthonous individuals-as evidenced in our study-can provide ample diversity in allele combinations, upon which natural selection can act. Our results provide important insights into how wild populations may genetically benefit from human-mediated translocations that will undoubtedly become increasingly necessary as more and more species face elevated extinction risks due to human activity and global climatic changes.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available on GenBank at https://www.ncbi.nlm.nih.gov/genba nk/, accession numbers MW646490-MW646518 (mitogenomes) and BioProject number PRJNA702749 (transcriptome raw data).