Ischnura damselflies were collected during the summer of 2015, 2016, 2019 or 2020 depending on their geographical origins. Individuals were caught in the field and stored in 95% ethanol in a -20°C freezer until further analysis. Specimens include 87 individuals from seven local populations in South Sweden (15 males and 72 females of all three morphs), and 105 other individuals from twelve other geographic regions, including Finland (Åland islands and mainland), Scotland, France, Cyprus, Greece, Ukraine, Belgium, the Netherlands, Italy, Montenegro and Japan (Table 1). The seven Swedish local populations are all located within a few kilometres from each other; thus, samples were grouped under a unique ‘Southwest’ population for the rest of the study. Similarly, samples from Belgium and the Netherlands were grouped as one unique population, denoted as ‘Vinne-Walem’. In contrast, the samples from isolated parts of Italy were separated in three populations based on their relative geographical locations, and denoted either ‘Northern’, ‘Central’, or ‘Southern’ (Figure 1, Table 1).
We also included 20 specimens belonging to three other Ischnura species in this study:
- Thirteen specimens of pumilio, collected in 2019 from four Swedish populations. Ischnura pumilio co-occurs sympatrically with I. elegans in this region but falls in another major phylogenetic clade of the Ischnura tree, and is less closely related to I. elegans than are the following two species .
- Three specimens of genei, an allopatric species to I. elegans endemic of the western Mediterranean region. Two samples of I. genei were collected from the two insular populations of Sardinia and Sicily (Italy), and the third sample is from Corsica (France) 
- One unique specimen of saharensis collected from Morocco. The species is also an allopatric species to I. elegans and is distributed across North Africa.
Specimens were dissected in sterile conditions to avoid cross specimen contamination. We extracted the DNA from the abdomen of each damselfly, except for the Japanese samples, for which DNA was extracted from one leg, following the protocol of a Qiagen DNeasy Blood & Tissue Extraction Kit (Qiagen, USA). The quality of all DNA extracts was tested by PCR, through the amplification of the 5’-end region (~654bp) of the cytochrome oxidase I (COI) mitochondrial gene using the primers LCO-1490/HCO-2198 designed by Folmer et al. . Only samples that were positive for the COI amplification were included in the following analyses.
All sequences were deposited into the GenBank databases (Accession #on-going submissions). In total, we amplified and sequenced four mitochondrial regions (COI, COIb, COIIa and NDI) to test for a selective sweep, one nuclear marker (PRMT) to test for population bottlenecks, and two Wolbachia genes (ftsZ and wsp) to characterize strain diversity. Note: an extra Wolbachia gene (fbpa) was sequenced for few wEle1-infected specimens from Sweden (primers details in Table S2). Purified PCR products were sent to Macrogen (Macrogen Europe, Inc.) for single strand direct forward Sanger sequencing. All sequences were manually curated and aligned using Geneious Prime 2020.2.4 (https://www.geneious.com) and AliView . Double peaks in the chromatograms were treated as either evidence of contamination (for mtDNA), multiple infections (for Wolbachia DNA), polymorphism (for mtDNA and nuclear DNA), or sequencing noises (all). Wolbachia and mtDNA sequences showing such patterns were not included in the following analyses, while the analysis of those double peaks in the nuclear locus sequences allowed us to identify heterozygotic and homozygotic specimens at polymorphic sites (see below).
Ischnura elegans nuclear haplotype diversity
To identify diversity at the nuclear level, we isolated 400bp of the successfully sequenced nuclear gene PRMT from 48 specimens (27 uninfected and 21 Wolbachia-infected) from 10 populations ([AB], [RK], [IK], [SW], [VC], [MI], [NI], [SI], [CY] and [GR]), with 2 to 12 specimens per population. However, Italy ([MI], [NI] & [SI]) and Japan ([RK] & [IK]) carry Wolbachia strains that are divergent from wEle1, and which may have altered the genetic of those populations in ways that are impossible to fully test with the current dataset. Therefore, we did not include these sequences in further analyses. The final sample size was 31 sequences, including 23 Wolbachia-uninfected specimens and 8 Wolbachia-infected specimens from five populations (Table 1). To test whether the nuclear gene from Wolbachia-infected and uninfected damselflies has evolved under neutrality, we performed neutrality tests in DnaSP v6.0  by calculating Tajima’s D (Tajima, 1989) and Fu and Li’s F  metrics, and by estimating nucleotide diversity (π) and haplotype diversity (Hd) (Table 2). We also estimated the observed heterozygosity levels at two nuclear polymorphic sites for both the Wolbachia-infected and uninfected specimens.
Wolbachia and mitochondrial haplotype diversity, phylogenies and haplotype networks
All Wolbachia Sanger-sequences were BLAST-ed against the Wolbachia PubMLST database (https://pubmlst.org/wolbachia/)  to find the corresponding or closest alleles at each locus. Additionally, all five MLST genes (ftsz, gatB, hcpA, coxA, and fbpA) and the wsp gene of the strain wEle1 were also extracted from the whole genome project of a Swedish I. elegans . All Wolbachia reads were identified and isolated from the raw read data of the I. elegans genome project . The wEle1 assembly was built by mapping reads to two previously sequenced Wolbachia genomes (wPip  and wMel ) using bwa mem version 0.7.8 . The properly mapped pairs were extracted using samtools 1.8. The isolated Wolbachia paired reads were assembled into a draft genome using spades version 3.9.0 at kmers 21,33,55,77, and 99. The wEle1 draft genome assembly (Nscaffold= 893; N50= 5 523bp; longest scaffold=53 331bp, genome size= 1.4MB) is available as supplementary fasta file.
We concatenated sequences of the ftsZ and wsp genes from each Wolbachia strain characterized in this study, and from five additional strains (wMel, wRi, wClec, wBm, wPip; Isolate id number: 1, 11, 36, 37, 1808 in Wolbachia pubMLST database, respectively) previously assigned to the A-, A-, F-, D-, and B-Wolbachia supergroups, respectively. The phylogenetic analyses were conducted in IQ-Tree on XSEDE  implemented in CIPRES , using the genes separately (Figure S1) or concatenated (Figure 2). ‘Model Selection’  was selected to allow for the search of the best model in CIPRES. The partition type was set to allow the two partitions (one for each gene) to have different speeds . The best fit substitution models were decided by running ‘-m TESTNEW’ in IQ-Tree. Bootstrapping was conducted using ‘Ultrafast’ and ‘SH-aLRT’ bootstrap methods (Hoang et al. 2018) in IQ-Tree with 1000 replicates. The ‘TN+F+I’ and ‘TPM3+F+G4’ models were applied to ftsZ and wsp genes, respectively, as the best fit models with highest BIC (Bayesian information criterion) scores. All other setting options were left as default. The pairwise genetic distances between Wolbachia strains were calculated in MEGA-X . The best phylogenetic trees were visualized in FigTree (http://tree.bio.ed.ac.uk/software/figtree/), and rooted using the wBm-D and wClec-F strains as outgroups (Figure 2).
Additionally, we built two types of mitochondrial haplotype networks: (A) one based on the COI 5’-end region only (598bp), and (B) a second based on all four mitochondrial regions (2375bp). The networks were built using POPART  with the median joining method  (Figure 3). To the mitochondrial sequences produced by the present study, we added mitochondrial sequences from the same markers from any species of the I. elegans clade (i.e. I. elegans, I. genei, I. saharensis, I. graellsii and I. fountaineae), publicly available in GenBank before July 2020 (Table S3). Note: only sequences with a length equal or longer to 600bp were included to ensure the performance of the analyses.
Wolbachia selective sweep
The estimation of the timing of the Wolbachia sweep in I. elegans was carried out following the method described by Rich et al. . We first estimated the neutral mutation rates at the third position of fourfold and twofold synonymous codons of the four mitochondrial genes separately. The open reading frames of mitochondrial genes were found by blasting the nucleotide sequence against the mitochondrial proteome of I. elegans . The mitochondrial sequences of I. elegans and I. pumilio were aligned in MEGA X  in order to calculate the number of nucleotide differences and the number of fourfold and twofold synonymous sites between the two species. Jukes-Cantor correction  was applied to correct for multiple substitutions. Lastly, the neutral mutation rates were calculated based on the divergence time between I. elegans and I. pumilio, estimated between 10.4 to 21.7 My before present  (calculations were repeated twice, using each extreme of that divergence time range). Consequently, the age of the infection can be estimated using the following equation:
Our estimation was only based on mitochondrial haplotypes that were carried by the infected individuals. With this method, the corrected number of substitutions was estimated as 137.91 = 245(-3/4)ln[1 – (4/3) x (97/245)] and 111.11 = 332(-1/2)ln[1 – 2 x (81/332)] among fourfold and twofold degenerate codons, respectively. Assuming the maximum estimate divergence time at 21.7 Mya , we estimated that neutral mutation rates on fourfold and twofold synonymous sites would be expected to be 1.30% and 0.77% per site per million years, respectively. If the minimum estimate of 10.4 Mya was assumed, the neutral mutation rates are estimated to be 2.71% and 1.61%, respectively. These estimations are biologically reasonable if we assume the general divergence rate of mtDNA in arthropods at 1.1-1.2% per site per million years (Brower 1994).
Finally, to test whether the mitochondrial genes from Wolbachia infected and uninfected damselflies have evolved under neutrality, we performed neutrality tests in DnaSP v6.0  by calculating two population genetic statistics: Tajima’s D (Tajima, 1989) and Fu and Li’s F , and by estimating nucleotide diversity (π) and haplotype diversity (Hd).
Ischnura elegans colour polymorphism, sex differences, and Wolbachia infection in Sweden
All statistical analyses were performed in R version 3.6.1 . We used Chi-square test to investigate the association between Wolbachia infection and sex, and between infection and colour morph in female I. elegans. As most populations had a limited and incomplete sampling per sex and colour morph (e.g. only one female from Åland and from Finland mainland; only two morphs present in the Japan, Cyprus and Scotland populations, see Zenodo open data submission doi: 10.5281/zenodo.4445061), this test was only performed using the Swedish specimens. Fisher’s exact test was applied as an improvement of Chi-square test when the expected value of any cells of the contingency table is below five (Table 3).