Quantitative trait locus analysis of parasitoid counteradaptation to symbiont-conferred resistance
Host and parasitoid linesBlack bean aphids (A. fabae) were reared on their host plant Vicia faba (Fabaceae) in a climate chamber at 22 °C with a 16-h photoperiod to ensure clonal reproduction. Two sublines of A. fabae clone A06-407 were used: the original A06-407 clone, which was free of any known defensive endosymbionts, and the modified A06-407 clone harboring the H. defensa strain H76 (Vorburger et al. 2009; Dennis et al. 2017). The original (H. defensa negative) and the modified (H. defensa positive, harboring the H. defensa strain H76) aphid lines are in the following called H− and H+, respectively.We used two experimentally evolved populations of the parasitoid wasp L. fabarum. One was adapted to the presence of Hamiltonella in host aphids, the other was not. Wasp populations were established by Dennis et al. (2017) from a mixture of nine collections of sexually reproducing, haplo-diploid L. fabarum from six locations across Switzerland. Experimental evolution was conducted by rearing wasps exclusively on H− or H+ aphids, leading to counteradaptation in the H+ treatment; wasps reared on H+ aphids evolved an improved ability to parasitize H+ aphids compared to wasps reared on H− aphids (see Dennis et al. 2017 for more details). After maintaining treatments for 24 generations in 4 replicate populations each, replicates were combined and treatments were continued unreplicated at a population size of 200 individuals (see Rossbacher and Vorburger 2020 for details). Until the onset of the experiments presented here, parasitoid populations had been reared for approximately 140 generations on either H− or H+ aphids (since September 2013). At this point, the population reared on H+ aphids was able to parasitize H+ aphids nearly as well as H− aphids, whereas the population reared on H− aphids was only able to parasitize H− aphids but not H+ aphids. In the following, we refer to the wasp population adapted to H+ aphids as R ( = Resistant to Hamiltonella) and to the population adapted to H− aphids as S ( = Susceptible to Hamiltonella).Experiment 1: characterization of general inheritance patternsTo determine whether the evolved ability to parasitize H+ aphids is mainly determined by the larval or the maternal genotype, and whether it shows a dominant or recessive inheritance pattern, crossing experiments were combined with no-choice bioassays over two generations of wasps. In the first generation, all possible combinations of males and females from the R and S populations were crossed in order to quantify their ability to reproduce on H+ aphids (Table 1). Assuming that this ability is governed by a single Mendelian locus with two alleles (R and S), which are fixed in the respective populations (likely an oversimplification), allowed us to postulate three mutually exclusive hypotheses (H1–H3) that make different predictions for the outcome of these crosses (Table 1). To indicate genotypes and ploidy, crosses are depicted in the following as, e.g., RR × S, meaning that a (diploid) female from the R population was crossed with a (haploid) male from the S population.Table 1 Prediction of female offspring survival and reproduction in experimental crosses of evolved Lysiphlebus fabarum populations under three different hypotheses.Full size table(H1) The counteradaptation is larval and dominant. Under H1, RR × R, RR × S, and SS × R crosses are expected to produce female offspring on H+ aphids, as homozygous RR and heterozygous RS female larvae would be of the R phenotype and thus counteradapted. Homozygous SS daughters from SS × S crosses would fail to develop. If the counteradaptation was larval but inherited in an intermediate rather than dominant fashion, the expectation remains the same as under H1, albeit with the possibility that RR × S and SS × R crosses produce fewer female offspring than RR × R crosses.(H2) The counteradaptation is larval and recessive. Under H2, only RR × R crosses would produce female offspring on H+ aphids. RR × S, SS × R, and SS × S crosses are expected to not produce any female offspring as their heterozygous (RS) or homozygous (SS) daughters would be of the S phenotype and thus not counteradapted.(H3) The counteradaptation is maternal. Under H3, the RR × R and RR × S crosses are expected to produce female offspring and the SS × R and SS × S crosses are not, as the genotype of the mother is decisive for offspring survival. If both maternal and larval effects were at play, the sex ratio in offspring from the RR×S crosses is expected to be male biased compared to RR × R crosses, due to a disadvantage of RS larvae compared to RR larvae, while haploid male larvae have an R genotype in either case.To isolate wasps prior to use in experiments, mummies (parasitized aphids approaching parasitoid emergence) were collected and stored individually in 1.5 ml Eppendorf tubes. Thus, adult wasps had never encountered another wasp or aphid before (naive virgins). Zero-to-3 days after hatching, the wasps were paired and given 20–120 min for mating in 1.5 ml Eppendorf tubes. Although there was no control whether mating occurred in the given amount of time, mating was usually observed within the first 30 s of having wasp pairs in the same tube. Then the wasps were released on a caged plant with an aphid colony consisting of a known number of 0–48-h-old H+ aphid nymphs. The mean ± standard deviation (SD) number of aphid nymphs provided per cross was 43.5 ± 14.9. Adult wasps were removed from colonies 24 h after release. Nine days after adding wasps, plants were enclosed in cellophane bags and left to dry out at 22 °C for hatching and subsequent sexing and counting of wasp offspring. Differences in numbers of female offspring between the different crosses of the first generation were analyzed with Mann–Whitney U tests. A generalized linear model (GLM) was used to analyze differences in sex ratios. Statistical analyses were performed using R version 3.5.2 (R Core Team 2018).Because findings from the first generation of crosses supported H3 (see “Results”), two extensions of H3 (H3.1 and H3.2) were tested in a second generation of crosses to determine whether the maternal counteradaptation was dominant or recessive (Table 1). To this end, we tested the ability of 20 virgin female offspring from 10 RR × S crosses (i.e., heterozygous RS females) to reproduce on H+ aphids. The mean ± SD number of aphid nymphs provided per RS female was 21.9 ± 9.5.(H3.1) The counteradaptation is maternal and dominant. Under H3.1, RS females are expected to reproduce successfully on H+ aphids, because they are of the R phenotype. They are expected to produce only male offspring as they are virgins (arrhenotokous parthenogenesis). This scenario is indistinguishable from cytoplasmic inheritance, which would require further examination.(H3.2) The counteradaptation is maternal and recessive. Under H3.2, RS females are not expected to reproduce on H+ aphids, because they are of the S phenotype.Experiment 2: crosses and phenotyping for QTL studyTo obtain a mapping population and phenotype data, a crossing scheme similar to the one by Pannebakker et al. (2011) was realized (Fig. 1). The crossing design relied on two main assumptions: First, we assumed that the alleles responsible for the counteradaptation are fixed in alternative states in the R and S populations. Second, due to the findings from the first experiment, we assumed the counteradaptation to be recessive and determined by the maternal genotype (see “Results”). In the first generation (P generation), a single S female was crossed with an R male to produce heterozygous female RS offspring (F1 generation). F1 females were allowed to reproduce as naive virgins to produce a recombinant male-only mapping population (F2 generation, Fig. 1A). F2 males were then backcrossed into the R background (each male with one RR female) to produce F3 female offspring for phenotyping (Fig. 1B). All reproduction up to the emergence of F3 females took place on H− aphids (Fig. 1) to avoid any selection. P individuals, F1 females and F2 males were stored in 1.5 ml Eppendorf tubes at −80 °C for subsequent genotyping.Fig. 1: Experimental crossing procedure for QTL analysis.Crossing design used to obtain a F2 mapping population for genotyping (A) and a F3 population for phenotyping (B). In a first step, two P generation individuals (parents), a diploid female from the symbiont-susceptible population, and a haploid male from the symbiont-resistant population were crossed to obtain 17 heterozygous F1 hybrid females. F1 hybrid females were allowed to reproduce as virgins—i.e., arrhenotokous parthenogenesis—to obtain 354 recombinant F2 males (mapping population), which were either carrying the S (susceptible) or the R (resistant) genotype. Recombinant F2 males were backcrossed with females of the resistant population to produce semi-recombinant F3 females. Sister F3 females have identical chromosomes of paternal origin and are thus considered clonal sibships. Two hundred and forty-four clonal sibships consisting of one to two sister F3 females were allowed to reproduce as virgins on a colony of symbiont-protected (H+) aphid hosts for phenotyping. Bar colors represent genomic regions originating from different parental populations and letters under sex symbols indicate the ploidy levels and genotypes.Full size imagePhenotyping was conducted by letting naive virgin F3 females oviposit for 24 h on colonies with a known number of approximately 24–72-h-old H+ aphid nymphs and subsequently counting their offspring as previously described. The average ± SD number of aphid nymphs provided was 40.9 ± 13.6. Wasps were added to the aphid colonies in an open Eppendorf tube. If possible, two sister F3 females from the same recombinant F2 father were added to each aphid colony in order to reduce the occurrence of false negatives, i.e., random failures to reproduce that are unrelated to the females’ genotype, e.g., due to harmful handling or death before oviposition. F3 sister females are identical concerning their paternal chromosome set and share the same R population background concerning their maternal chromosome set. They are considered clonal sibships (Pannebakker et al. 2011).The phenotype we measured was the number of wasp offspring produced per H+ aphid colony. This measure exhibited strong variation and zero inflation. To improve its value as a proxy for counteradaptation, the measure was corrected for certain variables in the phenotyping set-up that could have influenced offspring production independent of the F3 genotype. We used the zeroinfl function of the R-package pscl (Zeileis et al. 2008) to fit the following full model by zero-inflated Poisson regression:n_offspring ~ n_nymphs + n_wasps_added + all_removed + any_found_dead + any_in_tube | n_nymphs + n_wasps_added + all_removed + any_found_dead + any_in_tubewhere n_offspring is the number of offspring wasps produced, n_nymphs is the number of aphid nymphs, i.e., potential hosts, provided, n_wasps_added is a factor describing whether one or two wasps were added to the aphid colony, all_removed is a factor describing whether all wasps could be recovered 24 h after adding them to the aphid colony, any_found_dead is a factor describing whether any of the wasps were dead after 24 h, and any_in_tube is a factor describing whether any of the wasps were found in the tube rather than on the plant after 24 h. Parameters before and after the | symbol are components of the Poisson and the zero-inflation part of the model, respectively. The full model was reduced to a minimal model by performing backwards elimination with the function be.zerofinl from the R-package mpath (Wang 2020). The final minimal model was:n_offspring ~ n_nymphs + all_removed + any_in_tube | n_wasps_added + any_in_tube.Residuals of the minimal model were used as the corrected count phenotype for QTL mapping. We also assessed offspring presence presence/absence as an additional binary phenotype. Due to its simplicity, a binary phenotype may be less prone to environmental variation and more appropriate if counteradaptation is a Mendelian trait.DNA extraction and sequencingDNA extraction from 354 F2 males, 17 F1 females, and the two P individuals was performed adapting the LGC-sbeadex Livestock D protocol (LGC Genomics, Berlin, Germany). In addition to these experimental individuals, 30 wasps from an asexual, isofemale line of L. fabarum (line CV17-84) were processed to quantify genotyping error. Due to their mode of reproduction and maintenance at small population size, CV17-84 individuals are expected to be genetically nearly identical. F2, F1, and P individuals and three pools of 10 CV17-84 wasps each were crushed in liquid nitrogen prior to lysis. Extraction from individual samples was downscaled and included the following adaptations: lysis was done with PN buffer during 2 h at 60 °C with 1:10 protease solution, the lysate was incubated with binding mix during 20 min and elution was done at 60 °C. Extraction from pooled samples was, besides doubling the amount of protease, done following the manual. DNA concentration of each sample was measured using a Spark 10 M Multimode Microplate Reader (Tecan, Switzerland). Quality of DNA obtained with the used protocols was tested on a Nanodrop spectrophotometer (Thermo Fisher Scientific, USA) and on agarose gels. ddRAD library preparation was adapted from the protocol by Peterson et al. (2012). Restriction enzymes MfeI and TaqI were used for double digestion of up to 50 ng DNA per sample. After ligation of barcoded adapters to each individual sample, samples were combined in 12 pools with 24–36 samples each. Eleven pools contained one sample of 50 ng DNA from the CV17-84 wasps and 23–35 other samples (F2, F1, or P). Fragment size selection was performed on each pool with AMPure XP beads (Beckman Coulter, USA) (0.6× and 0.09×) and followed by selection of biotinylated P2 adapters. This was followed by PCR with KAPA HiFi HotStart ReadyMix (Roche, Switzerland) to amplify DNA and add 12 different Illumina primers to identify pools. Pools were then purified and combined into a final library. Mean fragment size of the library was 606 bp, as measured with the 2200 TapeStation (Agilent, USA), which corresponds to a mean insert size of 470 bp. The library was sequenced in a single lane of an SP flow cell on an Illumina NovaSeq 6000 System with 2 × 150 bp paired-end sequencing (at Functional Genomic Center, Zürich). A total of 307.2 million paired-end reads were obtained from P, F1, and F2 individuals and 11 CV17-84 control samples.GenotypingWe used the dDocent pipeline (Puritz et al. 2014; Puritz et al. 2014) for genotyping. Reads were demultiplexed with the process_radtags function of the STACKS package (v 2.14, Catchen et al. 2013) with disabled filtering of degraded cut sites, which led to 304.2 million demultiplexed paired-end reads. BWA-MEM (v 0.7.17, Li and Durbin 2010) was used with default settings to map reads to the reference genome of L. fabarum (Lf_genome_V1.0.fa, Dennis et al. 2020). On average (±SD), 1.585 (±1.058) million reads were assigned per sample during demultiplexing. out of which an average of 81.66% were mapped and retained after filtering for mapping quality (Supplementary Table S1). We called 547,092 variants using freebayes (v 1.3.1, Garrison and Marth 2012) with the default settings from dDocent pipeline specifying population (corresponding to the generation P, F1, F2, or CV17-84) and ploidy of individuals. The VCF-file was then split into a dataset containing 355 haploid individuals, i.e., males (one P, 354 F2) and a dataset with 29 diploid individuals i.e., females (1 P, 11 CV17-84, 17 F1). The dataset with diploids was filtered following the dDocent filtering pipeline up until removing indels, retaining 2456 single-nucleotide polymorphisms (SNPs). The following changes were made to the tutorial: the minimum quality score (–minQ) was set to 20, the minimum mean depth (–min-meanDP) was set to 10, and the maximum mean depth (–max-meanDP) was set to 400. The haploid dataset was then transformed to allelic primitives and filtered to contain only the 2456 SNPs that were retained in the diploid dataset. The VCF files containing haploid and diploid samples were then transformed to SNP tables using samtools (v 1.9, Li et al. 2009) and custom bash scripts. A custom R-script was then used to filter the SNP tables and create an input file for linkage mapping with MSTmap (Wu et al. 2008). The retained SNPs are homozygous in the mother, biallelic among the two parent individuals, and known in both parent individuals. Additionally, we tested for segregation distortion, removing SNPs that deviate significantly from an allele frequency of 50% based on a chi-square test with Bonferroni-corrected false-discovery rate of 5%. For each allele in each offspring (F2) male, alleles were recoded as “A” for maternal, “B” for paternal, and “U” for unknown. SNPs missing in >50% of individuals and individuals with >50% unknown genotypes were removed. The dataset used for linkage mapping contained 351 F2 individuals and 1838 SNPs of which 3 were removed by MSTmap internal filters leading to a final dataset of 1835 SNPs contained in the linkage map.Quantification of genotyping errorGenotyping error rate was quantified by counting mismatches between the supposedly identical genotypes of 11 CV17-84 DNA samples that were sequenced as part of 11 different pools. The 1835 SNPs used for QTL mapping were used as a template to filter SNPs in the dataset with CV17-84 individuals with vcftools (–positions flag). A SNP table containing CV17-84 genotypes was then analyzed in R to quantify genotyping error. For each pair of CV17-84 samples, the proportion of genotype mismatches was counted and averaged over all comparisons to obtain an estimate of mean genotyping error. Unknown genotypes were not counted as mismatch. The mean percentage of pairwise mismatches among the 11 CV17-84 samples ranged from 0.8392 to 1.706% with an average of 1.207%. The average mismatch measure was employed as an estimate for the genotyping error during analyses with R/qtl (Broman et al. 2003).Linkage map and QTL mappingLinkage mapping was performed with MSTmap (Wu et al. 2008) using the following settings: population_type = DH, distance_function = kosambi, cut_off_p_value = 0.000001, no_map_dist = 15.0, no_map_size = 2, missing_threshold = 0.25, estimation_before_clustering = no, detect_bad_data = yes, objective_function = COUNT. The resulting distance matrix was processed with R to contain only marker locations, Linkage group (LG) ID, and map distance. The new linkage map was edited in order to use the same LG IDs and orientations as in the linkage map by Dennis et al. (2020).Phenotype data, genotype data, and the new linkage map were merged with a custom R script to produce an input file for R/qtl (Broman et al. 2003). After reading the dataset with R/qtl, its cross type was transformed to recombinant-inbred by selfing (convert2riself function) because this expects no heterozygotes and genotype frequencies at 0.5, which fits our crossing scheme. We tested for duplicated genotypes ( >90% similarity between individuals), checked for switched markers using the checkAlleles function, and plotted recombination fractions (Supplementary Fig. S1), none of which indicated any problems. Intermarker distance was estimated with the est.map function, setting map function to “kosambi” and tolerance to 10−4. The resulting map was used as new linkage map with cM as map unit. Conditional genotype probabilities were calculated at a step size of 0.1 cM. The scanone function was used to calculate logarithmic of the odds (LOD) scores over the genome using the default (EM) algorithm with nonparametric and binary model for the corrected count phenotype and the additional binary phenotype, respectively. Significance thresholds were calculated by conducting 1000 permutations and choosing a 5% cut-off corresponding to the significance threshold at an alpha of 5%. The 95% approximate Bayes confidence interval was then calculated for the chromosome with significant LOD score. After simulating genotypes 1000 times with a step size of 0.1 cM and pulling genotype probabilities at the peak LOD, the explained phenotypic variance was estimated with the fitqtl function.Candidate gene identificationAs RADseq loci are usually short and represent a small proportion of the genome, they are unlikely located in candidate genes themselves. The 95% approximate Bayes confidence interval of the single significant QTL we identified includes all markers on scaffold tig00000002, upwards of 311,170 (bp). Thus, we considered tig00000002 from position 311,170 on as region for searching candidate genes. Gene annotations were retrieved from the recently published L. fabarum genome (Dennis et al. 2020). In addition, we identified putative venom and toxin genes in the L. fabarum genome in order to explore this function among candidate genes. To do so, we collected venom protein sequences from several parasitoid wasp species: Nasonia vitripennis (Danneels et al. 2010), Chelonus inanitus (Vincent et al. 2010), Microplitis demolitor (Burke and Strand 2014), Fopius arisanus (Geib et al. 2017), Diachasma alloeum (Tvedte et al. 2019), Cotesia congregata (Gauthier et al. 2021), Leptopilina boulardi, Leptopilina heterotoma (Goecks et al. 2013), and Aphidius ervi (Colinet et al. 2014); and retrieved candidate animal toxin proteins (7151 sequences) from the UniProt Animal Toxin Annotation Program database (UATdb, Jungo et al. 2012). These proteins were then matched to L. fabarum proteins by blastp (-e-value More
