The effect of genetic differentiation on reproductive isolation is assessed by measuring gametic compatibility and hybrid fitness using three population replicates for both HL and BL.
Populations and rearing conditions
We used three human-associated populations (London, UK, collected and introduced to culture in 2008; Budapest, Hungary, 2010; Watamu, Kenya, 2010) and three bat-associated populations (Hanušovice, CZ, 2016; Hoštějn, CZ, 2016; Raškov, CZ, 2016). The human-associated population replicates were chosen with respect to the large distances of their places of origin. The choice of the bat-related populations was limited to those well habituated to the artificial feeding system, where only populations from Czech Republic were available. All populations were reared in an incubator at 27 °C (optimal temperature to maintain weekly interval for feeding, egg laying, molting etc., according to our experience, or e.g. 26,32,43), at 70% relative humidity with a daily cycle of 12L: 12D. The populations as well as the experimental females and their offspring were artificially fed on a blood source that was novel for both HL and BL. The novelty is assumed based on a lipidomic analysis of bat blood, human blood, human blood conserved in CPDA (citrate phosphate dextrose adenine, Faculty Hospital Bohunice, Brno), or sperm of bedbugs fed with either of the three blood types showed clearly distinct profiles  [Additional file 1: Figure S5]. The CPDA-conserved blood was fed using parafilm bags and artificial feeding system . All populations had been habituated to the feeding system for at least two years prior to the experiment.
Design of the crossing experiments
Bed bugs entering the experiments were virgin. This was achieved by separating fed 5th instars individually into 96-well microplate wells, letting them molt into adults. Three-week-old females from each population were individually single-mated with a three-week-old male. The male came from any of the other five populations, that is, three between-lineage and two within-lineage crosses (N = 369; for sample sizes in individual lineage crosses, see Table 1). Within-population crosses were not executed to avoid confounding effects of inbreeding . The crosses were conducted in four batches across 18 months, all with approximately equal numbers of population cross combinations. The Raškov population was not available for the first batch.
The adult females were fed twice and mated immediately after the second blood meal. Males were fed twice with the second blood meal administered a week before mating. In this way, we ensured full sperm vesicles and males’ eagerness to mate . To ensure that the amount of sperm injected was similar, mating was standardized by interrupting after 60 s after successful intromission .
After mating, the females were isolated in a vial equipped with filter paper for egg laying. Females were fed weekly. We recorded if the females fed successfully and counted the number of eggs every week. To measure female fertility, fertilized and unfertilized eggs were distinguished, and the onset of infertility was determined following Otti et al. : fertile eggs are taut and whitish, with visible red eye spots of the developing embryo. Unfertilized eggs normally collapse soon after being laid and are greyish. The onset of infertility was established as the time point when an unfertilized egg was laid for the second time, to allow for one accidental fertilization failure. The total number of fertilized eggs was used to investigate the fecundity of the within- and between-lineage crosses.
If a female stopped laying eggs for two weeks in a row, we placed it in a well of a ventilated 96-well microplate. We recorded its survival in weeks and, since interspecific mating can be harmful in bedbug species [50, 51], female lifespan was used as an additional measure of the male effect on females. In total, we analyzed 369 females in 30 population combinations.
In order to analyze offspring fitness-related traits, the survival and body size, we collected 10–12 fed fifth instar nymphs from each female and transferred them to individual wells of 96-well microplates. This way we aimed to yield at least three sons and three daughters per female. After eclosion to adulthood, we recorded their survival without access to food. After the females and offspring perished, we measured their pronotum width as a representative scale of body size . In total, we analyzed 2791 offspring of 305 females.
Statistical analyses were carried out using RStudio 1.4.1717 (R version 4.1.1, ), , with packages lme4 , lmerTest , and coxme . First, we performed a Fisher exact test to investigate whether the number of females that did not lay any eggs differed between the lineage crosses (19 females did not lay any eggs: 5 females BL x BL (female x male lineage), 7 females BL x HL, 5 females HL x BL, 2 females HL x HL). Then we tested whether laying no eggs and the number of females dying during the egg-laying period differed between the lineage crosses (51 females perished during egg-laying: 7 females BL x BL, 8 females BL x HL, 23 females HL x BL, 13 females HL x HL). We analysed the first ten weeks of egg-laying because by then almost 95% of the females had stopped laying fertilized eggs (347 out of 369).
To analyze the total number and the number of fertilized eggs, we fitted linear mixed effects models (LME) with lineage cross type (between- x within-lineage cross) and female lineage (bat- x human-associated) including their interaction term as fixed factors and population cross (female x male population) and batch as random effects. Because females occasionally failed to feed in every week, the number of feedings varied between females. Therefore, we fitted the total number of feedings as a covariate in all models. Both fertilized and total number of eggs were significantly positively related to the number of times a female fed over the egg-laying period (LME: fertilized eggs—F1,354.59 = 61.970, P < 0.0001; total eggs: F1,354.56 = 95.746, P < 0.0001) (see also the supplementary results in Additional file 1).
For fertilized eggs, we also investigated whether egg-laying patterns differed between lineages by fitting an LME with fertilized eggs laid in each week as a response variable. The week, lineage cross type, female lineage, and their interaction terms were fitted as fixed factors. Using a binary variable of whether a female fed in a particular week, we accounted for variation in feeding behavior among females. Finally, we fitted individual, population cross, and batch as random effects.
For the analysis of the proportion of unfertilized eggs, we used the cbind() function in R to combine the number of fertilized and unfertilized eggs as a response variable. With this response variable, we then fitted a generalized linear mixed-effects model (GLME) with binomial distribution. Further, the lineage cross type, female lineage and their interaction term were fitted as fixed factors and the population cross and batch as random effects. Overdispersion was investigated using the DHARMa package . If overdispersion was detected, we accounted for it using an object-level random effect.
The coxme  and multcomp  packages were used to analyze the onset of infertility, female survival, and the survival of offspring. The lineage cross type, female lineage, and their interaction were fitted as a fixed factor, and the population cross and batch were fitted as random effects. Again, we fitted the number of feedings over the egg-laying period as a covariate. For the survival analysis of offspring (F1 adults), we additionally fitted sex and its interaction with the lineage cross type and the female lineage as fixed factors, and the population cross and mating pair were fitted as random effects. Because the offspring size could affect survival, we fitted the pronotum width as a covariate.
We analyzed adult offspring size, i.e., pronotum width, of the different populations fitting LMEs with female (mother) pronotum width, lineage cross type, female lineage, and offspring sex including their interaction terms as fixed factors, and population cross and mating pair as random effects.
We tested the independence and inbreeding levels of populations used in the study using a set of 9 microsatellite markers  (for the primer and multiplexing details, see Additional file 2: Table S5; for the PCR protocol, see Additional file 2: Table S6). We sampled 24–56 individuals from each population and extracted their DNA using PCRBioRapid kit (PB10.24, PCRBiosystems, London, UK). The fragment analysis was carried out at the Genetics Facility at the University of Bayreuth, using the Fragment Analyzer 5200 (Agilent Technologies, Waldbronn, Germany). Alleles were scored using PROSize version 3.0 . For the microsatellite data, we checked for the presence of null alleles using Micro-checker2.2.3 . The among-population Fst, allele frequencies, heterozygosity indices and AMOVAs (for each locus separately, 1000 permutations) within each host lineage and both lineages together were calculated using Genalex .