We introduce a new pattern of population genetic structure in a host-parasite system that can arise after secondary contact of previously isolated populations. Due to different generation time and therefore different tempo of molecular evolution the host and parasite populations reach different degrees of genetic differentiation during their separation (e.g. in refugia). Consequently, during the secondary contact the host populations are able to re-establish a single panmictic population across the area of contact, while the parasite populations stop their dispersal at the secondary contact zone and create a narrow hybrid zone. From the host’s perspective, the parasite’s hybrid zone functions on a microevolutionary scale as a “parasite turnover zone”: while the hosts are passing from area A to area B, their parasites turn genetically from the area A genotypes to the area B genotypes. We demonstrate this novel pattern on a model composed of Apodemus mice and Polyplax lice by comparing maternally inherited markers (complete mitochondrial genomes, and complete genomes of vertically transmitted symbiont Legionella polyplacis) with SNPs derived from the louse genomic data. We discuss circumstances that may lead to this pattern and possible reasons why it has been overlooked in the studies on host parasite population genetics.

Apodemus flavicollis mice were captured across the HZ of Polyplax serrata-specific lineages in 2018 and 2019 in the north-west of Czech Republic, using wooden snap traps. Mice were searched for lice by visual checking and combing. Lice were removed and stored in 100% ethanol in the -20°C. Genomic DNA from whole louse specimens was individually extracted using the Qiagen QIAamp DNA Micro Kit (Qiagen, Valencia, CA, USA). Twenty six lice ((15 SE and 11 SW)), 1 to maximum 3 from each parasitized host, were selected for genomic re-sequencing.

Whole genome re-sequencing was performed to generate metagenomic data used to map the SNPs, and to assemble genomes of the L. polyplacis symbionts and mitochondrial minichromosomes. gDNA libraries for twenty six louse specimen were constructed for paired-end Illumina sequencing with NovaSeq6000 instrument. All samples were sequenced on one Illumina Novaseq lane producing on average 62.5 million 150 bp paired-end reads (PE) per sample.

The contigs corresponding to Legionella polyplacis were visually identified using ORF prediction done in the Geneious package (the prokaryotic gene arrangement could be readily recognized by the density of predicted ORFs). In most samples, the genome of L. polyplacis was assembled into a single complete contig. In three assemblies (DPH41,19JA1_SW, 46MAN_SW) the symbiont genome was highly fragmented or could not be found, despite the good assembly quality of the louse genome (these samples were removed from the L. polyplacis analysis). Complete L. polyplacis genomes were annotated in RAST (Aziz et al., 2008) and aligned using Mafft algorithm implemented in Geneious. In three samples (Ne125b_Kot_SW, 99b_Pro_SE, 98c_Pro_SE) the rRNA region did not assemble correctly and was only represented by short fragments. To extend these fragments into a full length, we used the program aTRAM 2.0 (Allen, LaFrance, Folk, Johnson, & Guralnick, 2018). Phylogenetic analysis of the resulting 530,063 bp long matrix was performed in Phyml (Guindon et al., 2010).

Due to the conserved noncoding region shared by all minichromosomes, the Spades based assembly was not able to separate the different minichromosomes reliably. We, therefore, took an alternative approach. Using a single gene from each minichromosome as a reference (preliminary assignment of the genes to individual minichromosomes was based on the GenBank data available for P. spinulosa; acc. nos. KF647762-KF647771) we mapped the reads and extended the sequences in the program aTRAM 2.0 (Allen at al., 2018). To obtain annotations of the resulting sequences, we combined two methods. The first method utilized the web based server Mitos (Bernt et al., 2013). Since this method missannotated or entirely missed some of the genes, we corrected the results by a  blast based approach, taking advantage of the available annotated mitochondrial genome from P. spinulosa. We combined assembled minichromosome sequences into a custom database of P. serrata and used the 37 genes of P. spinulosa as blast queries (we run discontinous megablast and tblastx, both with E-value set to 10). We then aligned the P. serrata minichromosome sequences together with P. asiatica and P. spinulosa using Mauve (Darling, Mau, Blattner, & Perna, 2004). These alignments were used as a background for combining and correcting the annotations.