Iva Popovic
email: iva.popovic@uq.net.au

Nicolas Bierne
email: nicolas.bierne@umontpellier.fr

Federico Gaiti

Miloš Tanurdžić
email: m.tanurdzic@uq.edu.au

Cynthia Riginos
email: c.riginos@uq.edu.au

Pre-introduction introgression

Keywords: Mytilus, parallel evolution, introgression, hybrid zone, invasive species, genetic differentiation

come

Introduction

Biological invasions are leading threats to global biodiversity and central components of global change [START_REF] Chown | Biological invasions, climate change and genomics[END_REF][START_REF] Simberloff | Biological invasions: much progress plus several controversies[END_REF]. The impacts of introduced species (including invasive species) on native communities can be diverse and complex, and subsequently difficult to predict [START_REF] Oduor | Evolutionary responses of native plant species to invasive plants: a review[END_REF].

Since the earliest studies on biological invasions [START_REF] Baker | Characteristics and modes of origin of weeds. Characteristics and modes of origin of weeds[END_REF][START_REF] Baker | The evolution of weeds[END_REF], the literature now contains numerous examples of shared ecological traits contributing to increased survival and competitive abilities among introduced species (e.g., rapid growth rates and high fecundity through asexual reproduction; [START_REF] Richardson | Invasion science: the roads travelled and the roads ahead. Fifty years of invasion ecology: The legacy of Charles Elton[END_REF]. There is also growing evidence that local adaptation is often a significant factor in the success of introduced populations [START_REF] Colautti | Contemporary evolution during invasion: evidence for differentiation, natural selection, and local adaptation[END_REF][START_REF] Dlugosch | Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions[END_REF][START_REF] Prentis | Adaptive evolution in invasive species[END_REF]. Indeed, repeated evolution of ecologically important traits related to life history strategies, reproduction, growth and competition have been documented for introduced populations in both plants and animals [START_REF] Colautti | Contemporary evolution during invasion: evidence for differentiation, natural selection, and local adaptation[END_REF][START_REF] Ferrero | Invasion genetics of the Bermuda buttercup (Oxalis pescaprae): complex intercontinental patterns of genetic diversity, polyploidy and heterostyly characterize both native and introduced populations[END_REF][START_REF] Lee | Evolutionary mechanisms of habitat invasions, using the copepod Eurytemora affinis as a model system[END_REF][START_REF] Phillips | Invasion and the evolution of speed in toads[END_REF][START_REF] Turner | Rapid evolution of an invasive weed[END_REF][START_REF] Vandepitte | Rapid genetic adaptation precedes the spread of an exotic plant species[END_REF]. We know relatively less, however, about how heritable genetic variation is redistributed or altered to yield adaptive responses following introduction [START_REF] Bock | What we still don't know about invasion genetics[END_REF]. Because anthropogenic introductions are recent, selection on standing genetic variation from the native range is expected to be the primary genetic mechanism of adaptation in introduced populations [START_REF] Prentis | Adaptive evolution in invasive species[END_REF], thus increasing the likelihood that parallel genetic changes may accompany multiple independent introductions of the same species.

Genetic and genomic studies examining parallel evolutionary change in natural populations have shown that adaptation can be highly constrained and repeatable [START_REF] Bolnick | Non) parallel evolution[END_REF][START_REF] Conte | The probability of genetic parallelism and convergence in natural populations[END_REF][START_REF] Elmer | Adaptation in the age of ecological genomics: insights from parallelism and convergence[END_REF]. Replicated genetic changes across distinct lineages at various genomic scales are now documented, ranging from single amino acid polymorphisms (e.g., [START_REF] Hoekstra | A single amino acid mutation contributes to adaptive beach mouse color pattern[END_REF] to individual genes (e.g., [START_REF] Colosimo | Widespread parallel evolution in sticklebacks by repeated fixation of ectodysplasin alleles[END_REF][START_REF] Stern | Is genetic evolution predictable?[END_REF] and large genomic regions (e.g., [START_REF] Mérot | Intercontinental karyotype-environment parallelism supports a role for a chromosomal inversion in local adaptation in a seaweed fly[END_REF][START_REF] Van Belleghem | Evolution at two time frames: Polymorphisms from an ancient singular divergence event fuel contemporary parallel evolution[END_REF] including at the chromosomal scale (e.g., [START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF].

Comparative evidence for parallel evolution has also been demonstrated among closely related populations [START_REF] Jones | The genomic basis of adaptive evolution in threespine sticklebacks[END_REF]T. L. Turner, Bourne, Von Wettberg, Hu, & Nuzhdin, 2010), suggesting that adaptation can be predictable even on timescales relevant for invasive species research [START_REF] Marques | Genomics of rapid incipient speciation in sympatric threespine stickleback[END_REF]. Most studies of genetic repeatability, however, are predicated on the assumptions that focal populations have evolved under similar selective pressures [START_REF] Arendt | Convergence and parallelism reconsidered: what have we learned about the genetics of adaptation?[END_REF][START_REF] Bolnick | Non) parallel evolution[END_REF][START_REF] Rennison | Shared patterns of genome-wide differentiation are more strongly predicted by geography than by ecology[END_REF][START_REF] Stern | The genetic causes of convergent evolution[END_REF]. Invasive species are likely to experience a variety of new ecological conditions under which they did not evolve [START_REF] Colautti | Contemporary evolution during invasion: evidence for differentiation, natural selection, and local adaptation[END_REF] and may occupy vastly diverse abiotic conditions and habitats in their introduced ranges [START_REF] Guzinski | Population genomics of the introduced and cultivated Pacific kelp Undaria pinnatifida: Marinas-not farms-drive regional connectivity and establishment in natural rocky reefs[END_REF][START_REF] Tepolt | Adaptation in marine invasion: a genetic perspective[END_REF]. Invasive species are also likely to encounter congeneric species and gene flow may ensue if isolating barriers to reproduction are semi-permeable [START_REF] Viard | Marine invasions enter the genomic era: three lessons from the past, and the way forward[END_REF][START_REF] Viard | Anthropogenic hybridization at sea: three evolutionary questions relevant to invasive species[END_REF].

Introduced species adaptation frequently occurs in the context of gene flow between previously allopatric taxa [START_REF] Bay | Parallel introgression and selection on introduced alleles in a native species[END_REF][START_REF] Fitzpatrick | Rapid spread of invasive genes into a threatened native species[END_REF]. Such examples of very recent secondary contact can provide new insights into speciation genetics, especially how different genomic elements are shaped by selection and gene flow [START_REF] Dlugosch | The devil is in the details: genetic variation in introduced populations and its contributions to invasion[END_REF][START_REF] Ravinet | Interpreting the genomic landscape of speciation: a road map for finding barriers to gene flow[END_REF]. For instance, introgression from differentiated genetic backgrounds, such as native species, can contribute locally adapted alleles or new beneficial allelic combinations to introduced lineages [START_REF] Bay | Genomic islands of divergence or opportunities for introgression[END_REF][START_REF] Huerta-Sánchez | Altitude adaptation in Tibetans caused by introgression of Denisovan-like DNA[END_REF][START_REF] Martin | Interpreting the genomic landscape of introgression[END_REF]. Conversely, gene flow is expected to break down linkage disequilibria between sets of already well-adapted alleles, potentially impeding opportunities for post-introduction adaptation [START_REF] Kirkpatrick | Chromosome inversions, local adaptation and speciation[END_REF][START_REF] Nachman | Recombination rate variation and speciation: theoretical predictions and empirical results from rabbits and mice[END_REF][START_REF] Sousa | Understanding the origin of species with genome-scale data: modelling gene flow[END_REF]. Disassociation of locally adapted alleles can be slowed down, however, if beneficial alleles are tightly linked in low recombination regions, leading to heterogeneous genomic differentiation between populations in the presence of gene flow [START_REF] Bierne | Pervasive selection or is it…? why are FST outliers sometimes so frequent?[END_REF][START_REF] Harrison | Heterogeneous genome divergence, differential introgression, and the origin and structure of hybrid zones[END_REF]. These low-recombining regions likely experience strong linked selection during periods of allopatry that can accelerate lineage sorting and promote 'islands of differentiation' between closely related lineages [START_REF] Burri | Interpreting differentiation landscapes in the light of long-term linked selection[END_REF][START_REF] Burri | Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers[END_REF][START_REF] Cruickshank | Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow[END_REF]. Lowrecombining regions are also the most likely to accumulate differences associated with reproductive barriers or local adaption loci that resist or delay introgression during secondary contacts [START_REF] Barton | The barrier to genetic exchange between hybridising populations[END_REF][START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF][START_REF] Ravinet | Interpreting the genomic landscape of speciation: a road map for finding barriers to gene flow[END_REF][START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF].

Islands of genomic differentiation are common in natural divergence-with-gene-flow situations and emerging results suggest they are also frequent when secondary contact is human-caused. Genomic islands in regions of low recombination have been implicated in maintaining interspecific differences through association with introgression barriers between native-introduced oysters [START_REF] Gagnaire | Analysis of genome-wide differentiation between native and introduced populations of the cupped oysters Crassostrea gigas and Crassostrea angulata[END_REF]. Similarly, empirical evidence for parallel differentiation between replicated anthropogenic contact zones suggests that genome-wide admixture can lead to repeatable divergence when secondary contact is recent [START_REF] Bay | Parallel introgression and selection on introduced alleles in a native species[END_REF]Simon, Arbiol, et al., 2020). More broadly, theory suggests that interactions between concurrent divergent selection and gene flow may further constrain where differentiation can occur, thus favouring genomic 'reuse' or repeated differentiation in genomic regions (i.e. co-adapted alleles) that are less likely to be broken up by recombination [START_REF] Samuk | Gene flow and selection interact to promote adaptive divergence in regions of low recombination[END_REF]. Such genomic constraints may therefore increase the probability that the same loci are involved in maintaining differentiation between populations [START_REF] Samuk | Gene flow and selection interact to promote adaptive divergence in regions of low recombination[END_REF], despite the likelihood of diverse selective pressure experienced by introduced populations in novel environments [START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF][START_REF] Tepolt | Adaptation in marine invasion: a genetic perspective[END_REF].

Successful marine invasions provide exceptional opportunities to gain a better understanding of how gene flow and selection shape differentiation in introduced environments. Processes limiting adaptation rates in marine environments appear to be less severe than theory would predict [START_REF] Dlugosch | Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions[END_REF][START_REF] Nei | The bottleneck effect and genetic variability in populations[END_REF]. In general, high fecundity and highly dispersive larvae promote successive introductions of large numbers of individuals through human-mediated vectors (e.g., ballast water discharge; [START_REF] Carlton | Pattern, process, and prediction in marine invasion ecology[END_REF]. High propagule pressure is therefore a key mechanism suspected to help founding populations avoid the impacts of genetic drift and associated reductions in genetic diversity [START_REF] Blackburn | The influence of numbers on invasion success[END_REF][START_REF] Lockwood | The role of propagule pressure in explaining species invasions[END_REF][START_REF] Narum | Genomic patterns of diversity and divergence of two introduced salmonid species in Patagonia, South America[END_REF][START_REF] Roman | Paradox lost: genetic diversity and the success of aquatic invasions[END_REF]. Empirical studies in marine invasive species have supported these predictions, rarely demonstrating strong genetic bottlenecks in introduced populations [START_REF] Bernardi | Genomic signatures of rapid adaptive evolution in the bluespotted cornetfish, a Mediterranean Lessepsian invader[END_REF][START_REF] Gagnaire | Analysis of genome-wide differentiation between native and introduced populations of the cupped oysters Crassostrea gigas and Crassostrea angulata[END_REF][START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF][START_REF] Rius | Marine invasion genetics: from spatio-temporal patterns to evolutionary outcomes[END_REF], such that establishment could be aided by ample standing genetic variation available to selection [START_REF] Viard | Marine invasions enter the genomic era: three lessons from the past, and the way forward[END_REF]. Furthermore, the combination of frequent and dense introductions that likely underpin successful marine invasions ensure multiple worldwide introductions as (un)naturally replicated experiments of independent introduced populations occupying diverse coastal habitats alongside recent secondary contacts. Despite these advantages, few studies have investigated genome-wide variation and differentiation in marine non-native species across both native and introduced ranges (e.g., [START_REF] Bernardi | Genomic signatures of rapid adaptive evolution in the bluespotted cornetfish, a Mediterranean Lessepsian invader[END_REF][START_REF] Gagnaire | Analysis of genome-wide differentiation between native and introduced populations of the cupped oysters Crassostrea gigas and Crassostrea angulata[END_REF][START_REF] Guzinski | Population genomics of the introduced and cultivated Pacific kelp Undaria pinnatifida: Marinas-not farms-drive regional connectivity and establishment in natural rocky reefs[END_REF][START_REF] Riquet | Contrasting patterns of genome-wide polymorphism in the native and invasive range of the marine mollusc Crepidula fornicata[END_REF][START_REF] Rohfritsch | Population genomics shed light on the demographic and adaptive histories of European invasion in the Pacific oyster, Crassostrea gigas[END_REF][START_REF] Saarman | Introgression between invasive and native blue mussels (genus Mytilus) in the central California hybrid zone[END_REF][START_REF] Tepolt | Transcriptome sequencing reveals both neutral and adaptive genome dynamics in a marine invader[END_REF]. Inadequate resolution of population differentiation (e.g., [START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF] or complex demographic introduction histories (e.g., [START_REF] Rohfritsch | Population genomics shed light on the demographic and adaptive histories of European invasion in the Pacific oyster, Crassostrea gigas[END_REF], however, have precluded a detailed examination of whether introduced populations show parallel genetic responses to colonisation. Additionally, genomic studies of marine non-native species have either largely focused on genome scans of random genetic loci [START_REF] Bernardi | Genomic signatures of rapid adaptive evolution in the bluespotted cornetfish, a Mediterranean Lessepsian invader[END_REF][START_REF] Guzinski | Population genomics of the introduced and cultivated Pacific kelp Undaria pinnatifida: Marinas-not farms-drive regional connectivity and establishment in natural rocky reefs[END_REF][START_REF] Riquet | Contrasting patterns of genome-wide polymorphism in the native and invasive range of the marine mollusc Crepidula fornicata[END_REF][START_REF] Saarman | Introgression between invasive and native blue mussels (genus Mytilus) in the central California hybrid zone[END_REF]; but see [START_REF] Tepolt | Transcriptome sequencing reveals both neutral and adaptive genome dynamics in a marine invader[END_REF], or have not discriminated among putative sources of functional variation (i.e. protein-coding and non-coding variation) driving differentiation between native and introduced populations [START_REF] Tepolt | Transcriptome sequencing reveals both neutral and adaptive genome dynamics in a marine invader[END_REF]. Thus, we still have little insight regarding the genomic components of potentially selectively favoured variation that may underlie differentiation in introduced environments.

Here, we capitalise on a uniquely replicated series of introductions of the invasive marine mussel Mytilus galloprovincialis to investigate the role of recent evolutionary changes in the absence and presence of gene flow with a native congener. Within the native range of M. galloprovincialis, introgression with its sister species, Mytilus edulis, in Atlantic Europe, has led to pronounced genetic differentiation [START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF][START_REF] Quesada | Differential introgression of mitochondrial DNA across species boundaries within the marine mussel genus Mytilus[END_REF] and partial reproductive isolation [START_REF] El Ayari | The hidden side of a major marine biogeographic boundary: A wide mosaic hybrid zone at the Atlantic-Mediterranean divide reveals the complex interaction between natural and genetic barriers in mussels[END_REF] between Atlantic and Mediterranean M. galloprovincialis lineages. In the present study, we validate that each of these native lineages has been the source of multiple introductions in the northern and southern hemispheres (Fig. 1a). Moreover, there has been additional replication in the nature of receiving communities with respect to hybridisation: Atlantic M. galloprovincialis introduced into South Africa is not sympatric with any native congeners [START_REF] Branch | Can we predict the effects of alien species? A case-history of the invasion of South Africa by Mytilus galloprovincialis (Lamarck)[END_REF][START_REF] Daguin | Genetic relationships of Mytilus galloprovincialis Lamarck populations worldwide: evidence from nuclear-DNA markers[END_REF][START_REF] Grant | Mytilus galloprovincialis Lmk. in southern Africa[END_REF], whereas Atlantic M. galloprovincialis introduced into Australia hybridises with an endemic Australian taxon, Mytilus planulatus (Batemans Bay [START_REF] Popovic | Twin introductions by independent invader mussel lineages are both associated with recent admixture with a native congener in Australia[END_REF]. Similarly, Mediterranean M. galloprovincialis has displaced native Mytilus trossulus in southern California [START_REF] Geller | Decline of a native mussel masked by sibling species invasion[END_REF][START_REF] Hilbish | Historical changes in the distributions of invasive and endemic marine invertebrates are contrary to global warming predictions: the effects of decadal climate oscillations[END_REF][START_REF] Rawson | Hybridization between the blue mussels Mytilus galloprovincialis and M. trossulus along the Pacific coast of North America: evidence for limited introgression[END_REF] and only few hybrids are detected [START_REF] Saarman | Introgression between invasive and native blue mussels (genus Mytilus) in the central California hybrid zone[END_REF], whereas Mediterranean M. galloprovincialis introduced to Australia is admixed with M. planulatus (Sydney Harbour: Popovic et al., 2020). The unintentional introduction history of M. galloprovincialis provides an exceptional framework for investigating the reuse of standing genetic variation in the context of adaptation to new environments with contrasting opportunities for gene flow with native species.

We use transcriptome-wide markers and focus on repeated patterns of genomic differentiation to investigate whether introduced M. galloprovincialis in Australia and (non-hybridising) introductions in California and South Africa have undergone parallel genetic change relative to their native genetic backgrounds. We also quantify the extent of repeatability between independent admixture events in Australia involving divergent M. galloprovincialis lineages that differ in colonisation histories.

Comparisons of allele frequency deviations from the average genomic admixture rate allowed us to estimate genome-wide introgression patterns in each hybrid population and to determine whether loci differentiating parental taxa deviate towards introduced M. galloprovincialis ancestry or resist introgression (due to isolating barriers or selective processes). Because resolving the sources of selectively favoured variation can improve our understanding of the genetic mechanisms underlying adaptation [START_REF] Prentis | Adaptive evolution in invasive species[END_REF], we also examine whether differentiation and parallel responses involve both protein-coding transcripts and major non-coding elements of the transcriptome, long non-coding RNAs (lncRNAs). Non-coding elements of the transcriptome have been implicated as important sources of genetic variation during invasion [START_REF] Stapley | Transposable elements as agents of rapid adaptation may explain the genetic paradox of invasive species[END_REF] and there is evidence that higher RNA structures are involved with thermal adaptation in bacteria and may have similar roles in eukaryotes (de la [START_REF] De La Fuente | ncRNAs and thermoregulation: a view in prokaryotes and eukaryotes[END_REF][START_REF] Gaiti | Origin and evolution of the metazoan non-coding regulatory genome[END_REF][START_REF] Somero | RNA thermosensors: how might animals exploit their regulatory potential[END_REF]. Thus, lncRNAs may experience selection on ecological timescales relevant for studying populations or congeners [START_REF] Pang | Rapid evolution of noncoding RNAs: lack of conservation does not mean lack of function[END_REF][START_REF] Somero | RNA thermosensors: how might animals exploit their regulatory potential[END_REF], but remain unexamined in the context of biological invasions.

Our comparative investigation shows that repeatability of post-introduction differentiation and admixture varies between genetically distinct introduced lineages, where historical introgression influences contemporary evolution in independent introductions. The genetic data presented here represent the first genomic investigation of multiple M. galloprovincialis introductions at a global scale and the first comparative investigation of lncRNA evolution in an invasive species to date.

Methods

Sample collection, RNA extraction and sequencing

We obtained published RNAseq data for three Mytilus species: M. galloprovincialis from its native (n=15) and introduced range (n=18) in Australia, M. planulatus (n=5) and M. edulis (n=3) from our previous investigation of M. galloprovincialis introductions in the southern hemisphere (Table 1; Popovic et al., 2020;BioProject ID: PRJNA560413). We combined these data with new samples from introduced M. galloprovincialis in California (n=5) and South Africa (n=5) (Fig. 1a). Mussels were collected from rocky intertidal or subtidal environments (Table 1). Individuals in California were sampled outside of the southern hybrid zone boundary with Mytilus trossulus, to avoid sampling introgressed individuals [START_REF] Saarman | Introgression between invasive and native blue mussels (genus Mytilus) in the central California hybrid zone[END_REF]. All individuals were genotyped for the species diagnostic marker Glu-5' (Rawson et al., 1996) to confirm species identity. RNA was extracted from 10-20 mg of mantle tissue (preserved in RNAlater), using the RNeasy Plant Mini Kit (Qiagen, MD, USA) and following the Animal Tissues protocol. Strand specific cDNA libraries were constructed using the TruSeq stranded mRNA kit (Illumina) with average insert sizes of 250-300 bp. Paired-end libraries were sequenced across three lanes of an Illumina Hiseq2000 or across one lane of an Illumina Hiseq4000.

RNAseq data processing

We removed adapters from paired reads using Trimmomatic (v0.36) [START_REF] Bolger | Trimmomatic: a flexible trimmer for Illumina sequence data[END_REF].

We selected high quality reads using a Phred-scale average quality score of 20 within 4 bp sliding windows and a minimum size filter of 50 bp. Redundancy among high-coverage reads was reduced by digitally normalising each dataset using Trinity's insilico_read_normalization.pl script with a default kmer size of 25 and maximum read coverage of 50 [START_REF] Grabherr | Full-length transcriptome assembly from RNA-Seq data without a reference genome[END_REF]. Overlapping paired reads were merged using FLASH v1.2.11 [START_REF] Magoč | FLASH: fast length adjustment of short reads to improve genome assemblies[END_REF] with a minimum overlap length of 10 bp. Read mapping, lncRNA identification and all downstream analyses were conducted against a de novo M. galloprovincialis reference transcriptome obtained from a previously published study (Popovic et al., 2020 for de novo assembly details).

Identification of lncRNA transcripts

Putative lncRNA transcripts were identified following a published bioinformatic pipeline for lncRNA discovery [START_REF] Gaiti | Dynamic and widespread lncRNA expression in a sponge and the origin of animal complexity[END_REF]. Transcripts are classified as either protein-coding or non-protein coding using a series of computational tools to perform stepwise filtering, annotation and isolation of transcripts satisfying lncRNA criteria based on (i) homology to known proteins, (i) the presence of signal peptides, (iii) transcript length, (iv) open reading frame length, and (v) non-coding validation to retain the most likely non-coding candidate transcripts as predicted lncRNAs (details in Supplementary Material).

Mapping and identification of SNPs

We performed additional filtering of the reference assembly prior to population genomic analyses:

Transcripts showing high sequence similarity were clustered using Cd-Hit-Est [START_REF] Fu | CD-HIT: accelerated for clustering the next-generation sequencing data[END_REF][START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences[END_REF] with a minimum sequence identity of 95% of the shortest sequence. We removed contigs without a significant blastn hit to the M. galloprovincialis draft genome (e-value 10 -4 ; [START_REF] Murgarella | A first insight into the genome of the filter-feeder mussel Mytilus galloprovincialis[END_REF]. We also removed transcripts with significant blastn matches (e-value 10 -3 ) to the M. galloprovincialis male (Genbank reference: FJ890850.1, AY363687.2) and female (Genbank reference: FJ890849.1, AY497292.2) protein-coding and complete mitochondrial genomes. The resulting 'total' assembly of 143,093 nuclear sequences was used as a reference for Single Nucleotide Polymorphism (SNP) discovery (for analyses of parallelism) among four introduced and three nativerange M. galloprovincialis populations, M. planulatus, and M. edulis (Table 1). For comparisons of lncRNA contributions to population structure, we removed all remaining lncRNA loci from this reference assembly, leaving 106,333 transcripts considered as the 'full' transcriptome assembly. RNAseq reads were mapped to each reference assembly (i.e. total, full and lncRNA) using Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF]. PCR duplicates were removed using Picard MarkDuplicates (http://picard.sourceforage.net) and indexed BAM files were created with SAMtools. SNPs were called using Freebayes (https://github.com/ekg/freebayes), and variants were filtered with VCFtools [START_REF] Danecek | The variant call format and VCFtools[END_REF]. We removed sites below a genotype quality of 30 and a minimum mean depth coverage of 10 reads, singletons (minor allele count=2 across all populations) and indel variants, and all variants with missing genotype data. Genotypes were statistically phased using beagle v4.1 [START_REF] Browning | Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering[END_REF].

Analyses of population differentiation: Full transcriptome and lncRNAs

We visualised variation among individual genotypes using Principal Components Analysis (PCA) on SNPs derived from the full and lncRNA assemblies. Because we previously determined that introduced M. galloprovincialis in Sydney Harbour and Batemans Bay are introgressed with M. planulatus, we included M. planulatus samples from Tasmania as a reference parental population, but excluded M.

edulis. The full assembly resulted in 16,906 biallelic SNPs and we identified 471 biallelic SNPs for the lncRNA dataset following filtering. We calculated pairwise FST [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] for each dataset using the R package hierfstat [START_REF] Goudet | Hierfstat, a package for R to compute and test hierarchical F-statistics[END_REF] and tested for significance using the boot.ppfst function (10,000 bootstrap replicates). Given low differentiation between the eastern and western

Mediterranean populations (FST=0), we combined these populations for subsequent analyses.

To test the hypothesis that lncRNAs have a role in adaptation following introduction, we examined whether putative lncRNAs contribute to elevated differentiation among native and introduced populations relative to the full transcriptome. Alternatively, if lncRNA transcripts experience evolutionary constraints, we expected that transcripts would contribute significantly less to population structure compared to other expressed transcripts. To assess whether FST values were significantly different between full and lncRNA datasets, we randomly resampled 471 variants (i.e. number equivalent to lncRNA SNPs) from the full dataset 1000 times and calculated pairwise FST for each iteration to obtain a null distribution of estimates using custom scripts in R. We assessed whether the observed lncRNA FST values fell within the 5 th or 95 th quantiles of the null distribution to determine whether the lncRNAs contributed significantly to elevated or reduced population differentiation, correcting p-values for multiple comparisons at qvalue≤0.05 ('qvalue' R package; [START_REF] Dabney | qvalue: Q-value estimation for false discovery rate control[END_REF]. To explore diversity signatures that may indicate selection acting on lncRNAs, we examined the density of polymorphic sites (proportion of segregating sites per contig averaged across contigs) and contrasted observed and expected heterozygosity and allelic richness (using the smallest population size for rarefaction) for each dataset in hierfstat. PopGenome [START_REF] Pfeifer | PopGenome: an efficient Swiss army knife for population genomic analyses in R[END_REF] allowed us to calculate per locus within-population nucleotide diversity and Tajima's D values [START_REF] Tajima | Statistical method for testing the neutral mutation hypothesis by DNA polymorphism[END_REF] to provide an indirect indication of the number of rare variants. In all subsequent analyses of genetic parallelism (subsections below), we asked whether lncRNAs were among loci showing atypical genetic differentiation following introduction.

Patterns of genetic differentiation: Replicated introduced populations

To test the hypothesis that M. galloprovincialis introductions have undergone parallel differentiation relative to their native genetic backgrounds, we identified highly differentiated loci between introduced and native populations: (i) Batemans Bay and South Africa against Atlantic M. galloprovincialis and (ii)

California and Sydney Harbour against Mediterranean M. galloprovincialis. We calculated pairwise FST for each locus [START_REF] Nei | Molecular evolutionary genetics: Columbia university press[END_REF] using the function basic.stats() in hierfstat and converted negative FST values to zero. To focus analyses on long-term adaptive variation within individual contigs, we considered only SNPs that were shared between populations and also yielded the maximum FST value (maxFST) for each contig. If more than one shared SNP had a maxFST value, we retained a single SNP per contig to limit pseudoreplication of closely linked variants (e.g., [START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF]. We also repeated parallelism analyses using the maxFST value for each contig, irrespective of which SNP yielded the maximum value. Because FST measures are contingent on within-population diversity and are influenced by variation in minor allele frequency and sampled population sizes [START_REF] Berner | Allele frequency difference AFD-an intuitive alternative to FST for quantifying genetic population differentiation[END_REF], FST may not always reflect true allele frequency differences and may have decreased sensitivity at low differentiation levels, such as those characteristic of Mytilus populations; we therefore repeated analyses using absolute allele frequency differences between two samples (deltaP=|p1-p2|) to avoid assumptions associated with estimating maxFST.

If introduced populations diverge from their native genetic background in a repeatable manner or share the same history of colonisation (i.e. through stepping-stone migrations), we expected to see positive correlations in genetic differentiation between introduced-native populations, such that the most highly differentiated loci would overlap in different introduced locations. We also examined whether introduced populations originating from divergent genetic backgrounds in California (Mediterranean origin) and South Africa (Atlantic origin) displayed similar differentiation patterns from native reference populations.

For all joint comparisons, we used standardised major axis regression to estimate the strength of maxFST correlation. We estimated the slope and intercept (elevation) of the regressions using the SMATR R package [START_REF] Warton | smatr 3-an R package for estimation and inference about allometric lines[END_REF]. We categorised contigs as being highly differentiated in pairwise comparisons if they contained a maxFST value within the top 1% of the empirical distribution, whereby contigs within the 99 th percentile in both distributions were considered shared outlier contigs.

To explore whether our results are robust to sampling artefacts associated with small sample sizes and because some of our joint inferences involved pairwise comparisons of the same native population (e.g., Fig. 2a, Fig. 2b), we undertook a permutation-based test of significance to ascertain the effect of sampling error on the maxFST distribution and the influence of small samples size on outlier detection.

We took the total number of observed alternative allele counts per locus for each joint comparison and randomly resampled alleles within each locus and across all three populations (i.e. one native population and two introduced populations) of the same sizes as our empirical samples and recalculated FST and deltaP. This randomisation procedure (1000 permutations) allowed us to obtain the null distribution of summary statistics while (i) holding sample sizes and the number of loci sampled constant at the realised values and (ii) preserving the relationships across all three populations (i.e. one native population and two introduced populations). We recorded locus-specific estimates and used these values to generate marginal cumulative probability distributions to ascertain the probability of obtaining extreme values for some loci in joint comparisons due to random sampling effects.

Patterns of Introgression: Replicated contact zones

To evaluate genetic parallelism between independent admixture events in Australia, we quantified locusspecific introgression asymmetries as positive or negative deviations from the average genomic admixture rate, following the approach from Simon, Arbiol, et al. (2020). This allowed us to identify genomic regions either putatively resistant to introgression (positive deviations indicating an excess of introduced M. galloprovincialis ancestry) or permeable to introgression by heterospecific native alleles (negative deviations indicating M. planulatus ancestry). For each admixed population, we calculated the expected allele frequency of an allele (F exp ) as a function of the average population ancestry for the focal admixed population and its frequency in each parental reference population;

F exp = (f _Introduced * Q _Introduced ) + (f _Resident * Q _Resident );
where f is the frequency of the most common allele in the reference native-range population for the introduced species (i.e. f _Introduced in Atlantic or Mediterranean M. galloprovincialis) and its corresponding frequency in the resident species (i.e. f _Resident in native M. planulatus); and Q is the average admixture proportion in the admixed population corresponding to each putative parental ancestry component.

Given an expected allele frequency, we can calculate the locus-specific deviation from the average genomic admixture rate:

D = F obs -F exp ;
where F obs is the observed allele frequency in the admixed population.

We estimated genome-wide correlations in allelic deviations between independent contact zones to evaluate whether genome-wide departures from the genomic mean were similar between Batemans were included as predictors along with their interaction term.

For these analyses, population allele frequencies were calculated using VCFtools. We estimated individual ancestry proportions with the program ADMIXTURE [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF] to obtain average admixture proportions (Q) for each population. For ADMIXTURE analyses, we retained positions with up to 20% missing data across all populations, which resulted in 84,150 SNPs across 4941 contigs. We considered four putative ancestries corresponding to Atlantic and Mediterranean M. galloprovincialis and M. planulatus as a reference parental population, and M. edulis as an outgroup taxon. We ran ADMIXTURE with 100 iterations and used the cross-validation procedure with 50 replicates for K=4 genetic clusters. Because both Batemans Bay and Sydney Harbour showed >99% average ancestry belonging to M. planulatus and a single M. galloprovincialis lineage, we only considered two population ancestries for calculating expected allele frequencies.

It is well established that Atlantic M. galloprovincialis and M. edulis have a history of secondary contact, and that differential introgression with M. edulis contributes significantly to intraspecific outliers and partial reproductive isolation between Mediterranean and Atlantic M. galloprovincialis lineages [START_REF] El Ayari | The hidden side of a major marine biogeographic boundary: A wide mosaic hybrid zone at the Atlantic-Mediterranean divide reveals the complex interaction between natural and genetic barriers in mussels[END_REF][START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF]. Because divergent M. edulis genetic elements have been detected in Batemans Bay through secondary admixture with introduced Atlantic M. galloprovincialis (Popovic et al., 2020), we explored whether SNPs mapping to contigs strongly differentiating Mediterranean and Atlantic M. galloprovincialis lineages (contigs with maxFST ≥ 0.5) also showed elevated asymmetries in allelic deviations as the most likely candidates for reproductive barriers or local adaptation loci, where locus identity with respect to M. edulis introgression was treated as an additional variable in linear mixedeffect models (Boolean indicator for contigs with maxFST < 0.5 or ≥ 0.5).

Results

Predicted lncRNA transcripts

We predicted 44,096 putative lncRNA transcripts, representing 23.8% of the unfiltered reference de novo M. galloprovincialis transcriptome. This proportion is similar to high percentages reported for a marine sponge [START_REF] Gaiti | Dynamic and widespread lncRNA expression in a sponge and the origin of animal complexity[END_REF] and plants (garden pea, Kerr, Gaiti, Beveridge, & Tanurdzic, 2017; fern, Atallah, Vitek, Gaiti, Tanurdzic, & Banks, 2018) using the same computational pipeline, although this value is greater than the proportion of lncRNAs (14.6%) identified in the M. galloprovincialis digestive gland transcriptome [START_REF] Gerdol | RNA sequencing and de novo assembly of the digestive gland transcriptome in Mytilus galloprovincialis fed with toxinogenic and non-toxic strains of Alexandrium minutum[END_REF]. Despite stringent filtering, our final lncRNA assembly may include untranslated regions, other polyadenylated RNA transcripts that satisfied minimum length requirements (i.e. small regulatory RNAs, tRNAs and rRNAs) or novel Mytilus coding peptides with no similarity matches to protein databases.

Genetic differentiation: Full transcriptome and lncRNA datasets

Principal component analysis of 16,906 biallelic variants from the full assembly revealed genetic separation between introduced populations in California and South Africa (Fig. 1c), with the first and second PC axes explaining 6.51% and 5.69% variance among individual genotypes. This analysis also confirmed genetic structure previously identified between hybridising populations from Batemans Bay All populations displayed similar levels of observed heterozygosity, suggesting no reductions in genetic diversity compared to native populations that would indicate large bottleneck effects (Table S1, Table S2). S3). Compared to the full dataset, lncRNA variants returned lower FST values for most populations pairs (Table S3). Significant reductions in lncRNA FST, however, were only evident for some populations displaying high population structure in the full dataset (e.g., populations from divergent genetic backgrounds and comparisons with Tasmania; Table S3).

Considering all populations, we observed a lower density of polymorphic sites within lncRNAs (0.00397) compared to the full dataset (0.00453), although mean densities were not significantly different (Welch approximation t-test, P=0.07). There were no consistent differences between datasets in withinpopulation observed or expected heterozygosity and allelic richness that would suggest impacts of purifying selection on lncRNAs. Similarly, significant reductions in lncRNA polymorphism (Pi, withinpopulation nucleotide diversity averaged across loci) were only evident for the Atlantic population (P≤0.01, Table S1). lncRNAs did not show consistently lower Tajima's D values or diversity that would suggest stronger selective constraints or the possible action of balancing selection (i.e., positive Tajima's D; Table S4). Joint maxFST distributions revealed an effect of interspecific gene flow with M. planulatus on levels of differentiation between introduced and native populations. Admixed populations showed greater average maxFST against native-range populations from the Atlantic (Batemans Bay) (n=651 loci; standardised major axis regression, elevation=0.00721, P=0.00253, Fig. 2a) and Mediterranean (Sydney Harbour) (n=737, standardised major axis regression, elevation=0.00265, P=0.0115, Fig. 2b), compared to non-hybridising populations. Joint distributions were positively correlated in both comparisons. The maxFST correlation between Atlantic-derived populations was moderate, but significant (R 2 =0.532, P<0.0001; Fig. 2a), with highly differentiated contigs (FST>0.4) showing similar trends in both introduced populations (labeled in Fig. 2a). Three contigs were classified as shared outliers falling above the 99 th percentile of both empirical distributions and had joint maxFST values greater than expected by chance (P<0.001) based on random permutation tests that were conditional upon empirical samples sizes (Fig. S1). Correlations between Mediterranean-derived populations were also significant (R 2 =0.539, P<0.0001; Fig. 2b), with two contigs classified as shared outliers showing differentiation levels above the 99 th percentile in both introductions, although only one outlier was significant based on the simulated marginal P-value threshold (P<0.001) in both comparisons (Fig. 2b, Fig. S1). Joint distributions of deltaP estimates and maxFST per contig provided concordant results in all comparisons (Fig. S2, Fig. S3).

To investigate whether genetically divergent introductions may lead to similar genetic outcomes in the absence of hybridisation with native congeners, we tested for correlations in genetic differentiation between non-introgressed populations in California and South Africa against respective native reference populations (Fig. 2c). Joint maxFST comparisons resulted in fewer shared loci (n=318), a weak genomewide correlation (standardised major axis regression, elevation=0.0065, P=0.0043; R 2 =0.129, P<0.0001) and no shared outliers (Fig. 2c). In all comparisons, lncRNAs were poorly differentiated and none were differentiation outliers (Fig. 2). Incomplete and unordered scaffolds in the M. galloprovincialis draft genome assembly [START_REF] Murgarella | A first insight into the genome of the filter-feeder mussel Mytilus galloprovincialis[END_REF] precluded inferences regarding the genomic positions and the possibility of linkage between shared outlier contigs.

Repeatability of introgression patterns between contact zones

Mussels from Batemans Bay showed average ancestry proportions belonging to 35.2% Atlantic M.

galloprovincialis and 64.1% M. planulatus, whereas Sydney Harbour were composed of 73.3%

Mediterranean M. galloprovincialis and 26.7% M. planulatus (Fig. 1b). Joint comparisons of allelic deviations revealed heterogeneity in locus-specific introgression in both admixed populations, with some SNPs showing concordant positive or negative deviations; but there was no strong evidence for parallelism in genome-wide introgression patterns (R 2 =0.0003, P=0.02) (Fig. 3).

In Batemans Bay, parental FST (between M. planulatus and Atlantic M. galloprovincialis) predicted locus-specific allelic deviations in the direction of an excess of M. galloprovincialis ancestry (+D: slope=0.128, F 1,7421 =190.180, P<0.0001, Fig. 4a) with positive genome-wide deviations from the average admixture rate (intercept=0.0264, F 1,7421 =798.096, P<0.0001). Loci with maxFST ≥ 0.5 between Mediterranean and Atlantic populations had significantly higher deviation values (intercept=0.0463, F 1,7421 =4.002, P=0.0455; but no difference in slope, P=0.817) (Fig. 4a). Additionally, among three contigs (Contig15136, Contig34278, Contig33253) showing parallel genetic responses between Atlantic-derived populations (Fig. 2a), two are highly differentiated between parental Atlantic M. galloprovincialis and M.

planulatus and show negative deviations in Batemans Bay, pointing to possible selection against introduced alleles (Fig. 4a). Consistent M. edulis introgression contributing to intraspecific outliers between M. galloprovincialis lineages [START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF], SNPs with maxFST ≥ 0.5 between Mediterranean and Atlantic populations were largely fixed within M. edulis and show high allele frequency differentials (deltaP>0.6) for the M. edulis major allele (Fig. S4).

In contrast, allelic deviations in Sydney Harbour showed no effect of parental FST (between M.

planulatus and Mediterranean M. galloprovincialis) (slope=-0.0181, F 1,8678 =3.059, P=0.0803), although the mean genome-wide deviation was positive (intercept=0.0149, F 1,8678 =226.816; P<0.0001) (Fig. 4b).

Loci diagnostic of Mediterranean and Atlantic lineages showed significant negative deviations from the average admixture rate (intercept=-0.00344, F 1,8678 =4.76994, P=0.0290), but there was no difference in slope compared to all loci (P=0.928) (Fig. 4b).

In both admixed populations, the average deviation across SNPs localised to lncRNA transcripts was not significantly different than zero (Batemans Bay: n=15 SNPs, mean=0.332, P=0.148; Sydney Harbour, n=15 SNPs, mean=0.0124, P=0.432). Consistent with intraspecific comparisons (nativeintroduced M. galloprovincialis), lncRNA variants were generally poorly differentiated between M.

planulatus and native-range M. galloprovincialis (Fig. 4).

Discussion

Replicated species introductions provide opportunities to investigate repeated patterns of adaptation and local introgression across contact zones. In the present study, we clarify aspects of M.

galloprovincialis introduction history, validating evidence for an Atlantic origin for South African M.

galloprovincialis and Mediterranean-derived introductions in California [START_REF] Daguin | Genetic relationships of Mytilus galloprovincialis Lamarck populations worldwide: evidence from nuclear-DNA markers[END_REF] that complements our earlier results uncovering Atlantic and Mediterranean M. galloprovincialis admixture Australian mussels (Popovic et al., 2020). Select outlier loci showed repeated changes in allele frequencies from native reference populations, but there was limited evidence for strong genome-wide correlated shifts in allele frequencies between introduced populations. These results suggest that parallelism affects parts of the genome, with differentiation occurring predominantly within proteincoding regions, as lncRNA loci were minimally differentiated in pairwise comparisons. Comparing two independent hybrid populations, loci did not show parallel shifts in genome-wide introgression patterns from neutral expectations, implying different admixture dynamics when introduced lineages are genetically distinct. Interestingly, we observed that some loci with putatively M. edulis-derived alleles were especially resistant to introgression in Batemans Bay. This result is consistent with pre-introduction introgression shaping contemporary admixture whereby local interspecies introgression reinforces isolation between populations within species [START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF][START_REF] Foote | Killer whale genomes reveal a complex history of recurrent admixture and vicariance[END_REF][START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF]Simon, Fraïsse, et al., 2020, this issue). Together, our findings suggest that both selective processes and the demographic histories experienced by introduced lineages can influence the repeatability of post-introduction evolution.

Parallel differentiation among introduced populations at outlier loci

Replicated M. galloprovincialis introductions from the same native-range source lineage exhibited significant, albeit modest, genome-wide correlations (Fig. 2a, 2b), whereas introductions originating from genetically distinct lineages had much weaker correlations (Fig. 2c). Coincident with genome-wide evidence for parallel divergence, a small number of contigs were repeatedly involved in high levels of differentiation in Atlantic-derived populations (South Africa and Batemans Bay). Three contigs provided especially strong evidence for genetic parallelism following introduction, with high differentiation at the same nucleotide positions (maxFST>0.4) when considering joint outcomes across both population pairs in permutation tests (P<0.001) (Fig. 2a). Suggestively, Contig33253 returned a top blastn hit to the M. galloprovincialis hsp90-2 gene encoding heat shock protein 90 (HSP90); HSPs and other oxidative stress-related proteins have been shown to evolve under positive selection in M. galloprovincialis (Popovic & Riginos, 2020) and may also be targets of repeated selection within introduced environments or during long-distance transport experienced by migrant larvae in the ballast water of container ships.

Determining the evolutionary mechanisms by which parallel divergence arises across M.

galloprovincialis introductions is challenging, as a number of processes may converge towards similar differentiation patterns [START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF]. Parallel evolution is often interpreted as evidence for similar selection acting on the same genomic regions [START_REF] Bolnick | Non) parallel evolution[END_REF][START_REF] Butlin | Parallel evolution of local adaptation and reproductive isolation in the face of gene flow[END_REF], where standing genetic variation is adaptively reused [START_REF] Bock | What we still don't know about invasion genetics[END_REF][START_REF] Jones | The genomic basis of adaptive evolution in threespine sticklebacks[END_REF]. Genetic repeatability is also modulated by species-specific genetic architecture, either because complexes of adaptive alleles are protected in inversions [START_REF] Van Belleghem | Evolution at two time frames: Polymorphisms from an ancient singular divergence event fuel contemporary parallel evolution[END_REF] or because the differential effects of long-term linked-selection are conserved across similar recombination landscapes between species (e.g., Ficedula flycatchers, [START_REF] Burri | Interpreting differentiation landscapes in the light of long-term linked selection[END_REF][START_REF] Burri | Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers[END_REF]Helianthus sunflowers;[START_REF] Renaut | Shared selective pressure and local genomic landscape lead to repeatable patterns of genomic divergence in sunflowers[END_REF]Mimulus monkeyflowers, Stankowski et al., 2019;Heliconius butterflies, Martin et al., 2013). In the context of invasive species, background selection is less efficient at the earliest stages of differentiation (Burri, 2017) (i.e. between incipient native and introduced populations) and has limited power to generate differentiation in the absence of strong barriers to gene flow [START_REF] Stankowski | Widespread selection and gene flow shape the genomic landscape during a radiation of monkeyflowers[END_REF].

Consistent with this premise, genomic islands within low-recombining regions (those experiencing linked selection) also resist introgression between marine populations (e.g., European sea bass, [START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF]. Similarly, parallel divergence in the absence of ecological similarities, implicates an influence of recombination on both lineage sorting and introgression rate among other hybridising lineages (e.g., Pacific cupped oyster, [START_REF] Gagnaire | Analysis of genome-wide differentiation between native and introduced populations of the cupped oysters Crassostrea gigas and Crassostrea angulata[END_REF]; long snouted seahorse, [START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF]. Although genomic resources for M. galloprovincialis are insufficient to make inferences about whether loci diverging in parallel coincide with low-recombining regions, our findings suggest that similar selective pressures or genomic constraints linked to recombination landscapes or reproductive isolation may lead to repeatable differentiation in some parts of the genome. In addition, such forces can shape differentiation both in the absence or presence of interspecific gene flow with native taxa (Fig 2 ).

Differentiation estimates may be confounded by high allele frequency variance due to small sample sizes that can affect population genetic diversity and FST distributions [START_REF] Bierne | Pervasive selection or is it…? why are FST outliers sometimes so frequent?[END_REF][START_REF] Hoban | Finding the genomic basis of local adaptation: pitfalls, practical solutions, and future directions[END_REF]. While false positive outliers are unlikely to contribute to parallel differentiation by chance [START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF][START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF], permutation analyses allowed us to exclude the effects of small sample sizes as a substantial influence on the observed patterns, indicating highly significant (P<0.001) outliers, even when joint comparisons involved the same native population (Fig. 2a, Fig. 2b). Nevertheless, allele frequency shifts in native reference lineages may have contributed in part to correlated differentiation observed between introduced populations. Similarly, low genetic structure between native and introduced populations that is maintained by high standing genetic diversity and possibly ongoing propagule pressure, precluded our ability to rule out the possibility that shared colonisation histories (stepping-stone introduction routes) have led to parallelism at some loci. Indeed, difficulty in determining the native origins of introduced populations appears to be a ubiquitous issue in marine invasion studies and raises significant challenges for interpreting signals of differentiation associated with introductions [START_REF] Geller | Genetic perspectives on marine biological invasions[END_REF][START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF][START_REF] Viard | Marine invasions enter the genomic era: three lessons from the past, and the way forward[END_REF].

Despite largely positive correlations across M. galloprovincialis introductions, many loci showed location-specific shifts in allele frequencies, with little overlap in many differentiated contigs for both Atlantic and Mediterranean-derived joint comparisons (e.g., maxFST≥0.2; Fig. 2a, Fig. 2b). Introduced populations likely experience diverse selective pressures imposed by location-specific abiotic conditions, diverse community assemblages and the presence of semi-isolated native lineages [START_REF] Guzinski | Population genomics of the introduced and cultivated Pacific kelp Undaria pinnatifida: Marinas-not farms-drive regional connectivity and establishment in natural rocky reefs[END_REF][START_REF] Riquet | Contrasting patterns of genome-wide polymorphism in the native and invasive range of the marine mollusc Crepidula fornicata[END_REF]. In such cases, parallelism is expected to be reduced when different phenotypes are selected [START_REF] Bolnick | Non) parallel evolution[END_REF][START_REF] Thompson | Parallel genetic evolution and speciation from standing variation[END_REF].

Confidence in how location-specific selection shapes evolution in M. galloprovincialis introductions, however, will require additional empirical investigations.

Locus-specific introgression patterns differ between hybrid populations

Contrasting introgression patterns between hybrid populations in Australia revealed that gene flow can remodel population structure between introduced and native taxa in different ways depending on the introduced genetic background (i.e. Atlantic or Mediterranean source lineages). Specifically, joint comparisons of allele frequency deviations in Batemans Bay and Sydney Harbour indicated heterogeneity in locus-specific introgression rates, with no evidence for genome-wide correlations (Fig. 3), highlighting that introgression outcomes are not especially repeatable when divergent genetic lineages are involved. Locus-specific deviations were significantly correlated with parental FST in Batemans Bay (Fig. 4a), suggesting that the degree of divergence between hybridising species influences locus-specific introgression rates. Specifically, loci showed genome-wide deviations towards an excess of introduced M. galloprovincialis alleles introgressing into the M. planulatus native genetic background (Fig. 4a). In contrast, there was no evident relationship for Sydney Harbour mussels, although loci had a slight genomic tendency in the direction of M. galloprovincialis ancestry (Fig. 4b).

This result is unexpected as Sydney Harbour is a large and busy harbour likely subject to high propagule pressures of non-native larvae compared to the small docks and embayments of Batemans Bay, where we would expect M. planulatus to outnumber M. galloprovincialis (matching average ancestry proportions of 35.2% Atlantic M. galloprovincialis: 64.1% M. planulatus); thus, the predominant direction of introgression should be of native alleles into the introduced species background (i.e. negative deviation; [START_REF] Currat | The hidden side of invasions: massive introgression by local genes[END_REF].

Overall, our results support the hypothesis that differential introgression at some loci mirrors the divergence between historically allopatric taxa [START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF]. High rates of gene flow and recent divergence times typifying many marine hybrid zones may promote extensive introgression and divergence swamping between semi-isolated lineages, such that the history of divergence is only preserved at the most highly selected loci [START_REF] Bierne | Pervasive selection or is it…? why are FST outliers sometimes so frequent?[END_REF][START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF][START_REF] Riquet | Parallel pattern of differentiation at a genomic island shared between clinal and mosaic hybrid zones in a complex of cryptic seahorse lineages[END_REF].

In such cases, extrinsic selection may favour either invading or native alleles, depending on each hybrid zone environment. Alternatively, selection may be intrinsic, maintaining divergence in regions associated with reproductive incompatibilities. In Australian hybrid zones, a lack of concordance in admixture patterns among divergent introduced lineages is consistent with recent empirical work on M.

galloprovincialis introductions along the Atlantic coast of Europe (Simon, Arbiol, et al., 2020), where parallel patterns in allele frequency deviations were only observed when the same parental lineages were compared, whereas correlations were negative or non-significant when admixed populations comprised different M. galloprovincialis lineages. Thus, our results and those of Simon, Arbiol, et al. (2020), reveal that genomic backgrounds can strongly influence the outcomes of post-introduction admixture even among very closely related taxa.

Heterospecific alleles may drive genetic parallelism between Atlantic M. galloprovincialis introductions

From our analyses of paired introductions, we find that a small number of contigs were highly differentiated in both Atlantic-derived introduced populations (South Africa and Batemans Bay; Fig. 2a).

Among three contigs (Contig15136, Contig34278, Contig33253) showing the strongest evidence for genetic parallelism between Atlantic M. galloprovincialis-derived populations, two appear to resist introgression between Atlantic and Mediterranean M. galloprovincialis lineages (maxFST ≥ 0.5, contigs in green; Fig. 2a). Additionally, contigs with alleles putatively shared with Mytilus edulis (Fig. 5; Fig. S4)

were highly differentiated between parental Atlantic M. galloprovincialis and M. planulatus and showed negative deviations in Batemans Bay, pointing to possible selection against introduced genotypes (Fig. 4a). In Europe, local introgression from M. edulis into Atlantic M. galloprovincialis contributes significantly to differentiation at intraspecific outliers between Atlantic and Mediterranean M.

galloprovincialis [START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF]. M. edulis-derived alleles also appear to resist introgression between Mediterranean and Atlantic M. galloprovincialis lineages across a zone of secondary contact in Algeria [START_REF] El Ayari | The hidden side of a major marine biogeographic boundary: A wide mosaic hybrid zone at the Atlantic-Mediterranean divide reveals the complex interaction between natural and genetic barriers in mussels[END_REF] and can explain repeated patterns of intraspecific differentiation in some populations throughout the native range [START_REF] Bierne | Assortative fertilization and selection at larval stage in the mussels Mytilus edulis and M. galloprovincialis[END_REF][START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF][START_REF] Gosset | Differential introgression from a sister species explains high FST outlier loci within a mussel species[END_REF].

The influence of M. edulis introgression in naturally occurring native-range hybrid zones of Atlantic and Mediterranean M. galloprovincialis [START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF] is echoed in our results among introduced populations. In Batemans Bay, M. edulis alleles have undoubtedly accompanied introduced Atlantic M.

galloprovincialis (Popovic et al., 2020); when we restricted analyses of allelic deviations to variants within contigs strongly differentiating M. galloprovincialis lineages (contigs with maxFST ≥ 0.5, analogous to surveys using ancestry informative loci), parental FST had a significant effect on allelic deviations, with elevated mean departures from the genomic admixture rate (intercept=0.0463, P=0.0455) compared to analyses considering all loci (Fig. 4a). This pattern of asymmetric introgression illustrates an intriguing aspect of divergence between Atlantic M. galloprovincialis and M. planulatus, suggesting that highly divergent M. edulis alleles involved in reproductive isolation between M.

galloprovincialis lineages likely contribute to high locus-specific variation in introgression rates and may resist introgression with M. planulatus.

Consistently, we observed nearly fixed differences in two highly differentiated contigs in Atlantic-derived comparisons that are also putative M. edulis introgression outliers (Fig. 2a, Fig. 5). Contig15136

(encoding ribosomal protein L7a) was almost fixed for thymine (T) at position 965 in M. edulis (T) and in the Atlantic (T, except for one individual), while M. galloprovincialis and M. planulatus were fixed for cytosine (Fig. 5a). A similar pattern was observed in Contig34278, although the SNP frequency shift at site 872 is associated with one of M. edulis' alleles (Fig. 5b). The high maxFST values between Atlantic and South Africa (Fig. 2a) also likely involve hybridisation but of a different nature: mixed Atlantic and

Mediterranean ancestry proportions in one individual sampled in Cape Columbine, South Africa (Fig. 1), points to the possibility that introgressed genetic variation originating from the Mediterranean native range is likely to inflate maxFST, and could also reflect Mediterranean alleles resisting introgression with Atlantic M. galloprovincialis (that are heterogeneously introgressed by M. edulis), such that similar outcomes are evident in the presence (Batemans Bay) and absence (South Africa) of hybridisation with native congeners. For Sydney Harbor mussels, where partial M. edulis ancestry is not involved, parental divergence did not affect the departures of introgressing loci from neutral expectations (Fig. 4b).

From the circumstantial evidence that we have outlined, it appears likely that historical gene flow from M. edulis into Atlantic M. galloprovincialis (predating anthropogenic introductions) affects evolutionary dynamics in introduced mussel species, both in the context of hybridisation with the Australian M.

planulatus and the previously undocumented input of Mediterranean M. galloprovincialis to the invasive South African population (formerly thought to be exclusively Atlantic M. galloprovincialis-derived; [START_REF] Daguin | Genetic relationships of Mytilus galloprovincialis Lamarck populations worldwide: evidence from nuclear-DNA markers[END_REF]. Similar to the case of European port mussels, pre-introduction introgression with M.

edulis may have contributed to instances of repeatable admixture outcomes between introduced M.

galloprovincialis and the local genetic background (Simon, Arbiol, et al., 2020). However, given the small number of genotyped individuals in this study, we cannot make strong claims regarding introgression rates at individual loci. Sampling variance and genetic drift may lead to spurious allele frequency deviations, especially for loci where divergence between hybridising species is low. Likewise, genetic drift at individual loci may obscure parallel trends if only a few generations of backcrossing have occurred since the initial admixture event and there has been insufficient time for introgression to proceed (e.g., [START_REF] Gagnaire | Analysis of genome-wide differentiation between native and introduced populations of the cupped oysters Crassostrea gigas and Crassostrea angulata[END_REF]. Further investigations through simulations are warranted to better understand the dynamics of allelic deviation statistics under a spectrum of divergence scenarios and drift processes that may be important during invasive secondary contacts.

Based on genome-wide patterns, however, our findings in Batemans

Genetic differentiation: Long non-coding RNAs

Distinguishing among the major transcriptomic components of genetic variation underlying differentiation allowed us to characterise how lncRNAs shape population structure among native and introduced Mytilus populations (Fig 1c). LncRNA transcripts contributed significantly less to population genetic structure compared to variants derived from the full transcriptome (Fig. 1c; Table S1) and showed no evidence of high differentiation in analyses of parallelism (Fig. 2) or between congeners (Fig. 4). These findings suggest that lncRNA evolution may be constrained within lineages, with expressed protein-coding transcripts being more likely to dominate signatures of recent genetic differentiation.

Purifying selection on lncRNAs has been documented among Drosophila populations, and is consistent with suspected lncRNA functions in conserved developmental regulatory networks [START_REF] Gaiti | Long non-coding regulatory RNAs in sponges and insights into the origin of animal multicellularity[END_REF][START_REF] Perry | The functions of long noncoding RNAs in development and stem cells[END_REF]. However, measures of diversity among lncRNAs in M.

galloprovincialis indicated no consistent reductions that would support stronger purifying selection or balancing selection on these transcripts.

Despite notable progress in understanding the role of lncRNAs in speciation, the evolutionary processes influencing lncRNA transcripts are poorly understood, with mixed evidence supporting both sequence conservation and rapid evolution among lineages [START_REF] Gaiti | Long non-coding regulatory RNAs in sponges and insights into the origin of animal multicellularity[END_REF]. lncRNAs frequently evolve faster than protein coding sequences, as the primary sequence is not constrained by protein function [START_REF] Hezroni | Principles of long noncoding RNA evolution derived from direct comparison of transcriptomes in 17 species[END_REF][START_REF] Pang | Rapid evolution of noncoding RNAs: lack of conservation does not mean lack of function[END_REF][START_REF] Quinn | Unique features of long non-coding RNA biogenesis and function[END_REF][START_REF] Ulitsky | Evolution to the rescue: using comparative genomics to understand long non-coding RNAs[END_REF]; however, alternative modes of evolution, such as adaptive changes to secondary and tertiary structures may occur through nonnucleotide substitution mechanisms that cannot be detected by SNP-based measures of genomic differentiation as employed in this study. We acknowledge that lncRNA identification using the total assembly subjected to additional filtering (prior to population genomic analyses) resulted in a reduced number of lncRNA transcripts and polymorphisms. Thus, genotyping errors due to mapping biases or assembly errors in the absence of an annotated reference genome may have decreased detectable population structure and our results should be interpreted with some caution. For Mytilus mussels, we also cannot exclude the possibility that non-selective processes have contributed to patterns of lncRNA evolution and diversity observed in this study, including variation in mutation rate or evolutionary constraints linked to genomic synteny that do not operate in a population-specific way, or at least not within timescales relevant for biological invasions.

Conclusions

Here, we provide comparative evidence for parallel genetic changes in multiple M. galloprovincialis populations in the context of anthropogenic introductions and recent secondary contacts. Our study shows that species introductions, both in the presence and absence of gene flow with native congeners, may lead to similar genetic changes, especially for protein-coding loci. The extent of parallel differentiation and admixture patterns between contact zones, however, varied depending on the colonising M. galloprovincialis lineage. We observed stronger introgression asymmetries in contigs distinguishing Mediterranean and Atlantic M. galloprovincialis lineages compared to genome-wide trends in Batemans Bay, an observation consistent with M. edulis alleles influencing introgression barriers between introduced lineages (South Africa) and in secondary contact with native congeners (Batemans Bay). Genomic studies of marine invasive species and especially anthropogenic hybrid zones are still rare [START_REF] Viard | Marine invasions enter the genomic era: three lessons from the past, and the way forward[END_REF][START_REF] Viard | Anthropogenic hybridization at sea: three evolutionary questions relevant to invasive species[END_REF]. The lack of strong reproductive isolation among M. galloprovincialis and M. planulatus provides a contrast to secondary contact in Ciona tunicates [START_REF] Bouchemousse | Distinguishing contemporary hybridization from past introgression with postgenomic ancestry-informative SNPs in strongly differentiated Ciona species[END_REF][START_REF] Viard | Anthropogenic hybridization at sea: three evolutionary questions relevant to invasive species[END_REF] and oysters [START_REF] Gagnaire | Analysis of genome-wide differentiation between native and introduced populations of the cupped oysters Crassostrea gigas and Crassostrea angulata[END_REF] where reproductive isolation is strong. Thus, in introduced-native Mytilus hybrid zones there is the opportunity to observe the earliest stages of extensive admixture that may lead to consolidation into islands of divergence or to complete genomic swamping. Overall, our findings build on a growing number of investigations (e.g., [START_REF] Duranton | The contribution of ancient admixture to reproductive isolation between European sea bass lineages[END_REF][START_REF] Fraïsse | Local interspecies introgression is the main cause of extreme levels of intraspecific differentiation in mussels[END_REF]Simon, Arbiol, et al., 2020;[START_REF] Simon | How do species barriers decay? Concordance and local introgression in mosaic hybrid zones of mussels[END_REF] reporting significant impacts of historical introgression on the intrinsic genetic * Species identity and range for Mytilus samples is based on genome-wide analyses of species relationships (Popovic et al., 2020).

and

  Sydney Harbour, which showed intermediate placement between northern M. galloprovincialis populations and Australian M. planulatus (Fig. 1c; Popovic et al., 2020). ADMIXTURE analyses indicated shared ancestry between introduced M. galloprovincialis in California and Mediterranean populations; individuals sampled in South Africa clustered closely with Atlantic mussels, validating previous inferences regarding the origins of these introductions (Fig. 1b). PCA of genotypic variance confirmed that individuals sampled in California are not hybrids with Mytilus trossulus, where introgressed individuals would have shown divergent placement away from the Mediterranean cluster.

  lncRNA variant identification and filtering resulted in 200 transcripts with 471 SNPs. Lower expression among lncRNAs compared to coding transcripts (e.g.,[START_REF] Quinn | Unique features of long non-coding RNA biogenesis and function[END_REF] may have resulted in the removal of many variants due to low sequence coverage. Overall, PCA of lncRNA variants explained a lower proportion of variance across the first (5.92%) and second (5.10%) PC axes compared to the full transcriptome. For this dataset, only M. planulatus sampled in Tasmania and the admixed population in Batemans Bay showed genetic separation from all other populations which did not reflect clustering based on sampling region (Fig.1c). Pairwise comparisons of population differentiation for the full dataset indicated FST values significantly different than zero for all population pairs, with the lowest values for comparisons between the Mediterranean and California (FST=0.003; Table

Fig. 1 .

 1 Fig. 1. (a) Sampling locations for native (circle) and introduced (triangle) M. galloprovincialis populations and M. planulatus. Approximated distributions of introduced populations in Australia, California and South Africa are denoted by red (Mediterranean) and purple (Atlantic) indicating known native source lineages, and the M. planulatus (MP) range in Australia is shown in yellow. (b) ADMIXTURE results for K=4 genetic clusters, including M. edulis (ME). Each bar represents an individual belonging to one or more ancestral clusters, corresponding to different colours. (c) Principal components analysis (PCA) of M. galloprovincialis individuals sampled from geographic locations coloured by geographic location, where M. planulatus samples are included as a reference population. Top panel shows a PCA of 16,906 variants derived from the full transcriptome assembly and the bottom panel shows a PCA across 471 SNPs derived from lncRNA transcripts.

Fig. 2 .

 2 Fig. 2. Joint maxFST distributions for pairs of introduced M. galloprovincialis populations against native reference populations from the (a) Atlantic (elevation=0.00721, P=0.00253; R 2 =0.532, P<0.0001); and (b) the Mediterranean (elevation=0.00265, P=0.0115; R 2 =0.539, P<0.0001). (c) Joint maxFST distributions for introduced M. galloprovincialis populations from different native origins in South Africa and California (R 2 =0.129, P=<0.0001). For all plots, the grey dotted line corresponds to the 1:1 slope and the solid line denotes the major axis regression. The red dotted lines show the 99 th percentile for both distributions, above which contigs were considered outliers. Solid red lines indicate the probability (P<0.001) of observing outliers due to chance based on permutation tests. Green points correspond to contigs strongly differentiating Atlantic and Mediterranean M. galloprovincialis lineages (maxFST≥0.5) as putative loci with introgressed M. edulis alleles. None of these diagnostic loci were present in comparisons involving Mediterranean populations (Fig. 2b, 2c). Black open circles indicate lncRNA transcripts.

Fig. 3 .

 3 Fig. 3. Joint comparison of allelic frequency deviations from the average genomic admixture rate for hybrid populations in Australia. Mean genome-wide deviations are shown by the dotted lines for Batemans Bay (mean=0.00613) and Sydney Harbour (mean=-0.00320). Loci showing positive deviations from the genomic average indicate an excess of M. galloprovincialis ancestry from the

Fig. 4 .

 4 Fig. 4. Locus-specific allele frequency deviations for admixed populations as a function of parental divergence (snpFST) between M. planulatus (MP) and M. galloprovincialis (MG) lineages from the Atlantic and Mediterranean. (a) Significant effects of parental FST on allelic deviations in Batemans Bay (grey; slope=0.128, P<0.0001) suggest that differential introgression is influenced by the degree of parental divergence, with loci showing positive genome-wide deviations (intercept=0.0264, P<0.0001) towards an excess of M. galloprovincialis ancestry. SNP variants within contigs differentiating M.galloprovincialis lineages (green; contigs with maxFST ≥ 0.5) had significantly higher deviation values (intercept=0.0463, P=0.0455; but no difference in slope, P=0.817). (b) There was no effect of parental FST on allelic deviations in Sydney Harbour mussels (P=0.0803). The mean deviation for all loci was positive (grey; intercept=0.0149, P<0.0001), whereas diagnostic loci showed a slight negative departure from the average admixture rate (green; intercept=-0.00344, P=0.0290).

Fig. 5 .

 5 Fig. 5. Frequencies of polymorphic sites along (a) Contig15136 and (b) Contig34278. Graphs show allele frequencies in reference to the Mytilus edulis major allele. Focal SNPs with maxFST values for each contig and where signals of parallelism may be attributed to local introgression of M. edulis alleles are indicated.

  if highly differentiated loci deviate equally toward introduced or native ancestries (opposing deviations), we would not expect to see correlations with parental FST. For each admixed population, we compared allelic deviations against parental FST values (between M. planulatus and introduced M. galloprovincialis lineages) estimated for each SNP locus (snpFST). To account for nonindependence of loci on the same contig, we fit the data with linear mixed-effects models that included contig identity as a nested random effect using the maximum likelihood method and the lme() function in the R package nlme[START_REF] Pinheiro | nlme: Linear and nonlinear mixed effects models[END_REF]. Both parental FST and a Boolean predictor reflecting Mediterranean vs. Atlantic M. galloprovincialis differentiation (see below)

Bay and Sydney Harbour despite involving divergent M. galloprovincialis lineages. If admixture events have led to repeated genetic outcomes, we expected to see positive correlations between allelic deviations in separate localities. The repeatability of admixture can also be assessed by analysing introgression patterns independently within each admixture event. Specifically, if the degree of divergence between incipient species is inversely correlated to introgression rate

[START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF]

, we predicted that highly differentiated loci between native M. planulatus and introduced M. galloprovincialis lineages would display strong asymmetries in locus-specific introgression (either positive or negative deviations) relative to the average genomic admixture rate. Under the hypothesis of barriers to introgression or other sources of selection favouring either introduced or locally adapted alleles, genome-wide deviations and parental FST are expected to be correlated, whereby the most divergent loci show an excess of M. galloprovincialis ancestry (positive correlations) or tendencies towards resident ancestry (negative correlations). Under neutral expectations (no barriers to introgression) or

  Filtered datasets of native-introduced population pairs resulted in 11,906 (Sydney Harbour) and 10,088

	SNPs (California) for Mediterranean-derived populations, and 11,183 (Batemans Bay) and 8,902 SNPs
	(South Africa) for Atlantic-derived populations. Consistent with recent timescale of presumed M.
	galloprovincialis introductions into Australia (since the earliest records of European contact <600 years
	ago to 100 years ago with the advent of commercial shipping), maxFST distributions of shared loci
	between introduced and native populations were skewed towards low values, with few loci being highly
	differentiated between populations. Overall, introduced populations originating from Atlantic M.
	galloprovincialis were more differentiated from their native-range and showed greater genome-wide
	variance compared to Mediterranean-derived introduced populations. maxFST for individual contigs
	ranged between 0-0.6 (mean=0.0397, SD=0.0651, n=1634) for South Africa and 0-0.696
	(mean=0.0529, SD=0.0753; n=1779) for Batemans Bay. Introduced populations in California and
	Sydney Harbour were less differentiated from Mediterranean M. galloprovincialis, with maxFST ranging
	from 0-0.372 (mean=0.0325, SD=0.045, n=1729) in California and 0-0.375 (mean=0.0319, SD=0.0410,
	n=1873) for Sydney Harbour.

  Bay implicate partial reproductive isolation between Atlantic M. galloprovincialis and M. planulatus, which may also explain why we observed different outcomes when Mediterranean M. galloprovincialis was considered. This results suggests that different sets of incompatibility genes may be important for M. planulatus and divergent M. galloprovincialis lineages. While several reproductive isolating mechanisms are documented among several pairs of Mytilus congeners (i.e. asynchronous spawning, Gilg, Kirby, Sullivan, Knapp, & Hilbish,

	2007; differences in habitat preference, Gosling & McGrath, 1990, assortative fertilization, Bierne et al.,
	2002; Dobzhansky-Muller incompatibility substitutions at several loci, Bierne, Bonhomme, Boudry,
	Szulkin, & David, 2006), we currently have no information on the possible isolating mechanisms
	(ecological or genomic) between M. planulatus and other Mytilus species.

Table 1 .

 1 Summary of sampling locations.

	Taxon	Sampling Location	Range	Samples sequenced
				(RNAseq)
	Mytilus edulis	Darling Marine Station, Maine, USA	Native	3
	Mytilus galloprovincialis Primel, France (Atlantic)	Native	5
		Crique des Issambles, France (Western Mediterranean) Native	5
		Herceg Novi, Montenegro (Eastern Mediterranean)	Native	5
		Scripps Institute of Oceanography, California, USA	Introduced	5
		Cape Columbine, South Africa	Introduced	3
		Yzerfontein, South Africa	Introduced	2
		Sydney Harbour, New South Wales, Australia*	Introduced	9
		Batemans Bay, New South Wales, Australia*	Introduced	9
	Mytilus planulatus	Spring Bay, Tasmania, Australia*	Native	5

differences between marine lineages, which may also influence the repeatability of recent evolutionary change and the genetic outcomes of human-mediated admixture events.

Acknowledgements

We would like to thank F. Viard, J. Thia, D. Matias, M. James, and D. Ortiz-Barrientos for valuable discussions and insight on the manuscript. The QRIScloud computing cluster at the University of Queensland provided computational infrastructure for this project. This work is supported by the Australian Biological Resources Study (ABRS) National Taxonomy Research Grant (grant number RF216-11), and awards from the National Sciences and Engineering Research Council of Canada (to I.P), the Society for the Study of Evolution (to I.P) and the University of Queensland (to I.P).

Data Archive

RNAseq data is available on the NCBI Sequence Read Archive (SRA) (BioProject ID: PRJNA560413).

Raw RNAseq data for newly sequenced populations are deposited to the NCBI SRA (BioProject ID: PRJNA679634). Reference transcriptomes, lncRNA reference assemblies, and genomic datasets that support the findings of this study are openly available in the Dryad digital repository at https://doi.org/10.5061/dryad.xsj3tx9d9 and the Genomic Observatories MetaDatabase (https://geomedb.org).