For a schematic overview of the experimental design, see Fig. 2.
Experimental evolution
Protocol
The stock population (Stock population below) was allowed to expand for one generation and from this we established eight replicate experimental evolution populations, four selected for fighter morphs (F-lines) and four selected for scrambler morphs (S-lines). Each population was founded by 1,000 recently eclosed adults (500 random females and 500 random males of the desired morph). The classification of the morphs was based on visual inspection using a stereoscopic microscope and was unambiguous due to the discontinuous distribution of the phenotypes (Classifying male morphs below). Adults were allowed to interact freely for 6 days, all surviving adults (with previously laid eggs discarded) were transferred to a new container for 24 h of egg laying, after which adults were removed. The resulting offspring were allowed to mature over 13 days and 1,000 individuals from the newly eclosed adults selected for founding the following generation, again 500 random females and 500 random males of the desired morph, with this protocol repeated every generation (Extended Data Fig. 1). The isolation of nymphs to use virgins was unfeasible with our experimental design and population sizes. However, the period of 6 days after selecting the founders of the next generation and collecting eggs for the next generation was probably enough to displace most sperm stored by females mated with any unselected males due to the high number of remating that will be occurring over this duration (females on average remate after 80 min, ref. 88) and last male sperm precedence89. The timing of generation was chosen to reflect maturation rates from our stock population to avoid indirect selection on this trait. Moreover, a previous study90 showed there was no difference between male morphs in maturation rates and that over similar lengths of time to the protocol here the fertility of both morphs remains similar. Therefore, our protocol was not likely to impose strong differential selection on morph life histories.
Tracking morph proportion
We assayed the proportion of male morph in each population every 6–7 generations, by isolating 200 larvae (ten per vial) from the container, allowing maturation within vials and recording the morph of all males that eclosed (mean n = 86 per population, per generation, range 71–109). Our selection protocol was highly effective in driving an increase in the frequency of the desired male morph to >90% after 20 generations in both treatments, with this effect considerably faster within F-lines indicted by a significant two-way interaction between proportion of the desired morph and generation (χ2 = 39.9, d.f. = 6, P < 0.001; Extended Data Fig. 3). Selection for specific morphs was more effective than a previous experiment by Plesnar-Bielak et al.39, in which it took around 35 generations to drive both the desired morphs above 90%, we note that their lines selected for fighters also responded faster than those selected for scramblers. The difference between experiments in driving the frequency of desired morph >90% is probably a consequence of a longer interaction period (3 versus 6 days) in which the stored sperm of males before selection was able to be displaced and/or because selection was acting more efficiently in our larger populations. The difference between rate of changes in morph proportion between F- and S-lines in the current study, and also found by Plesnar-Bielak et al.39, may be associated with the genetic architecture of morph expression. Alternatively, selection could be less effective in scrambler lines if they are less efficient than fighters in displacing sperm of previous females’ partners, but this is unlikely as R. robini male morphs have previously been demonstrated to not differ in their sperm competiveness89.
Stock population
We established a stock population by mixing three laboratory populations that were collected from three sites in Poland (Krakόw, collected in 1998 and 2008, Kwiejce, collected in 2017 and Mosina, collected in 2017; Extended Data Fig. 1), where the line derived for material used in creating the reference genome (below) was also established from the same collections at Mosina in 2017. All populations were maintained in cultures with several hundred individuals per generation before mixing and establishment of the stock population. The mixing of distinct populations increased the genetic variance in the stock population, which otherwise would probably have been limited due to founder events and the limited population size of each of the contributing populations73, thus decreasing our power to detect the effects of SSTs on genetic variation. The newly mixed stock population was maintained with several hundred individuals per generation for roughly 12 generations before the onset of this experiment. This time period is probably enough to break linkage disequilibria that could have arisen due to mixing (for unlinked loci, linkage disequilibrium should decay by half each generation91).
One generation before establishing experimental evolution populations the proportion of male morphs was determined from 176 random males, indicating a roughly equal morph ratio (95 fighters, 81 scramblers) of the stock population (Extended Data Fig. 3).
General housing and husbandry
The stock population and experimental evolution populations were maintained in plastic containers (approximate length, 9.5 cm; width, 7 cm; height, 4.5 cm), filled with roughly 1 cm of plaster-of-Paris. The same containers were used when sampling mites for sequencing for the reference genome or resequencing from experimental evolution populations, but either replaced the plaster-of-Paris with 5% agarose gel or added a thin layer of 5% agarose gel above the plaster-of-Paris, respectively. The agarose gel was used to reduce the number of contaminates within our samples and on the basis of preliminary extractions that indicated that small pieces of plaster-of-Paris may reduce the quality of DNA during extractions. Individuals, pairs and small groups of ten mites were housed in glass vials (approximate height, 2 cm; diameter, 0.8 cm) and large groups of 60 or 150 mites in plastic containers (approximate height, 1.5 cm; diameter, 2 cm diameter or height, 1.5 cm; diameter, 3.5 cm diameter, respectively) all with an approximate 1 cm base of plaster-of-Paris. All plaster-of-Paris bases were completely soaked in water before mites were transferred into them. All mites were reared at a constant 23 °C, at high humidity (>90%) and were provided an excess of powdered yeast ad libitum.
Classifying male morphs
To illustrate the discontinuous distribution of the weapon and to demonstrate that this classification based on visual inspection is non-subjective, we performed phenotypic measurements from male mites from a population collected near Krakόw, Poland, that had previously been fixed onto microscope slides for a separate study66. The measurements taken were idiosoma (body without mouthparts) length and width of third proximal segment of the third right leg (genu). Measurements were preformed using Lecia DM5500B microscope and Lecia Application Suite v.4.6.1. We then performed an analysis to, first, determine whether the allometric relationship between idiosoma length and width of third pair of legs is best described as discontinuous and, second, to verify that classification by simple visual inspection matches the same classification from allometric analysis. One researcher performed all the measurements and classified each male as a fighter (n = 50) or scrambler (n = 50), a separate researcher was then given the measurements but not the classification of the male morph.
Broadly, guidelines for the analysis of non-linear allometries92 were followed. The log–log scatterplots of idiosoma length against leg width were visualized, which showed there was clear evidence for non-linear scaling relationships. Next histograms of idiosoma length, leg width and relative leg width (leg width/idiosoma length) were visualized (Extended Data Fig. 2a–c). Where a normal distribution of idiosoma length, and a binomial distribution in leg width and relative leg width are further indications of a discontinuous relationship. On the basis of the lowest point between the two peaks of the density plot of relative leg width (Extended Data Fig. 2c) males were classified as scramblers (relative leg width <0.125) or fighters (relative leg width >0.125). Replotting the log–log scatterplot of idiosoma length and leg width, and using the classification of morph described above clearly demonstrates the discontinuous allometric relationship of idiosoma length and leg width in R. robini (Extended Data Fig. 2d). Moreover, on the basis of the Akaike information criterion (AIC), the discontinuous model where males were assigned a morph (AIC = 646.5) clearly has a substantially better fit than a simple linear and quadratic models (AIC = 918.5 and 920.2, respectively). Further models were omitted from comparison (for example, breakpoint or sigmoidal) due to the clear discontinuous allometry observed. Finally, all 100 males were assigned the same morph by visual inspection and blind allometric analysis, demonstrating that the former is effective and accurate in classifying male morph.
Phenotypic assays
Fecundity assays were performed using experimental evolution females at F20 and F32. Eggs laid by females between days 4–8 were counted, encompassing the window of time of most evolutionary relevance for female fitness during maintenance of selection lines (that is, egg laying period in selection lines was between days 6–7) and also likely to capture variation in lifetime fecundity that remains largely consistent throughout the first 3 weeks of life93. Nymphs were individually isolated to gain virgin females, which on maturation females from each experimental evolution population (n = 30) were paired with a male from the stock population (15 with fighters and 15 with scramblers). Pairs were transferred to a new vial on day 4, with the pair being removed from the second vial after a further 4 days and all eggs in the second vial counted. If the male had died in the first vial, they were replaced with a stock male of the same morph. Any female deaths in the first or second vials were recorded.
Longevity assays were also performed at F20 and F32. At F20, females used in fecundity assays, including the stock male they were paired with (replaced if dead), were transferred to a new vial at day 8. After this point, vials were then checked every 2 days for female deaths and pairs were moved to new vials every 4 days. Males were replaced with stock males of the same morph if found dead. Similarly, at F20, on maturation males from experimental evolution populations (n = 30) were paired with stock females, vials were checked every 2 days and changed every 4 days, with females being replaced if dead. At F32, only female longevity was determined and was performed in groups; 30 experimental evolution females and 30 stock males (15 of each morph) were placed in plastic containers, two per experimental population. This logistically easier estimate of longevity was done due to local restrictions during the SARS-CoV-2 pandemic and the imposed limitations on people working closely together. Groups were checked for dead females every other day and all remaining live mites transferred to a new container every 4 days. When mites were transferred to a new container the sex and morph ratio were balanced to that of the remaining females, by either removing or adding males of the desired morph from the stock population.
To determine whether the survival of mites differed between F- and S-lines when competition between males was allowed, at F45 we created small colonies from each population and survival of males and females recorded over 6 days, the same period as used between selecting founders of the next generation and subsequent egg laying period. Colonies were at a 50:50 sex ratio, established with 150 newly eclosed mites placed into small plastic containers. This was approximately the same density after selection of the next generations founders during the maintenance of experimental evolution populations (150 mites in roughly 9.5 cm2 = 16 mites per 1 cm2; 1,000 mites in roughly 67 cm2 = 15 mites per 1 cm2). After 3 days, all colonies were checked and any dead mites identified by sex. After another 3 days, again dead mites were recorded and all surviving mites sexed and counted.
Additionally, at F45 we performed further fecundity assays to obtain estimates of inbreeding depression within experimental evolution populations. To establish family groups, larvae were isolated and on maturation F0 males and females (n = 16) from within the same experimental evolution population were paired together. Pairs were allowed to produce eggs for 48 h, after which adults were removed from vials. After hatching from each pair, 12 F1 larvae were isolated into new vials. On their maturation, these F1 mites were either paired with a full sibling, that is, from the same family, or with an individual from a different family but from the same experimental evolution population. When possible, we made two inbred and two outbred pairs with same family lines used. Again, pairs were allowed to produce eggs for 48 h before their removal for the vial. After a further 5 days, vials were checked for larvae, if larvae were present in the first vial six were individually isolated and the second vial discarded, if no larvae were present in the first vial the second vial was checked for larvae and, if present, they were isolated. This protocol therefore produced inbred and outbred individuals from within the same experimental evolution population. Which, as above, on maturation F2 inbred and outbred females were paired with stock males (fighter males only) and number of eggs laid between days 4 and 8 counted. Only a single female from each unique inbred or outbred family was used. Either due to pairs failing to produce offspring or there being no F2 females, samples sizes were not exactly equal. In total, 59 outbred and 55 inbred females from F-lines, and 56 outbred and 54 inbred females from S-lines were paired with stock males.
Phenotypic assay statistical analyses
All phenotypic analysis was conducted using R statistical software94 (v.3.5.2) and data were visualized using ggplot2 (ref. 95).
Analysis of male morph proportion was performed using a generalized linear mixed model with binomial error structure, fitted using lme4 (ref. 96). Where the proportion of desired morph was compared in model with morph selection and generation (as a factor) including their two-way interaction as explanatory variables, and population included as a random effect.
All fecundity data were analysed using generalized linear mixed models with Poisson error structures, fitted using lme4. Due to the differences in stock population males used between F45 and earlier generations, and slightly different rearing conditions between females in the fecundity assays from generations F20 and F32, they were analysed separately from data collected in F45. However, we noted that the fecundity of females in Fig. 5a was comparable to the outbred females in Fig. 5b. Explanatory variables fitted to fecundity data from F20 and F32 were, morph selection treatment, generation, including their two-way interaction term, and stock male morph. The explanatory variables fitted to fecundity data from inbreeding depression data were, morph selection treatment and status of female (that is, inbred or outbred), including their two-way interaction term. In both analyses, we included population as a random effect and an observation level random effect to account for overdispersion, we omitted fitting random slopes due to issues with increasing the complexity of random effects close to reaching a singular fit. Females that died before the end of the fecundity assay and those that laid zero eggs were removed from analysis. This excluded five females from F20 (three F-line and two S-line), 20 from F32 (13 F-line and seven S-line) and 16 from F45 (three inbred and three outbred F-line, and nine inbred and one outbred S-line).
Longevities of females at F20 and F32, and males at F20, were analysed separately using mixed effects Cox models, fitted using coxme97. In all analyses, we included a random effect of population, with morph selection treatment as an explanatory variable and extra variable of male morph included in female longevity analysis at F20. Survival of mites over 6 days at F45 was analysed using a GLM with counts of dead and surviving mites fitted with a quasibinomial error structure, the model included morph selection treatment and sex, including their interaction term, as explanatory variables. If individuals were lost due to handling error (that is, killed or escaped) they were right-censored during analysis.
Genome assembly
Sample origin
A line of R. robini originated from a wild-collected population from the Mosina region (Wielkopolska, Poland). In October 2017, onions were collected from the field and approximately 200 individuals of R. robini were identified under dissecting microscope. The line used for DNA isolation in the genome sequencing project was developed from full sib × sib mating for 14 generations (to maximize homozygosity) following and continuing the protocol described in ref. 67.
DNA extraction
For DNA extraction we used only mite eggs, that were laid by 500 females, collected in a container (see above for a description) Females were kept in this container for 3 days. After that time, they were removed, and eggs were filtered using fine sieves and washed for 1 min in 0.3% sodium hypochlorite solution and in Milli-Q water for 2 × 2 min to remove any potential foreign DNA contamination. These eggs were collected in 1.5 ml Eppendorf tube and after short centrifugation, the remains of the water (supernatant) removed with a pipette. The sample was immediately transferred to ice and prepared for DNA extraction. DNA was extracted using Bionano Prep Animal Tissue DNA Kit for HMW DNA isolation according to the manufacturer’s instructions. Briefly, eggs were smashed with a sterile pellet pestle on ice in 500 μl homogenization buffer; the sample was fixed with 500 μl cold ethanol and incubated 60 min on ice, after that time the sample was centrifuged at 1500g for 5 min at 4 °C and the supernatant was discarded. Next, after resuspension in a homogenization buffer pellet, this was cast in four agarose plugs as described in the original protocol. Agarose plugs were incubated with Proteinase K and Lysis buffer solution for 2 h with intermittent mixing. After that time, the digestion solution was replaced with a freshly made one and incubated overnight with intermittent mixing. According to the original protocol, after RNase A digestion and plug washing, DNA was recovered by incubation of the plugs in TE buffer, followed by plug melting and addition of agarase. Recovered DNA was dialysed and homogenized on a membrane for 45 min at room temperature and transferred to a clean tube with a wide bore tip.
Sequencing
Sequencing was done using Oxford Nanopore Technologies (ONT, MinION). Isolated DNA purified using AMPure XP beads and resuspended in H2O before library preparation. Two separate libraries were prepared using ligation sequencing kit, SQK-LSK109 and Rapid Sequencing Kit SQK-RAD004, respectively, according to the manufacturer’s protocols and were sequenced on a FLO-MIN106 R9.4.1 SpotON flow cell on a MinION Mk 1B sequencer (ONT). The total yields from sequencing were 484,700 reads (2,417,068,187 nt) with a read-N50 of 10,044 nt (ranging from 216,403 to 100). Base calling of the raw reads was done using Guppy (v.3.3) resulting in a total sum of the reads 7,979,616,172, equivalent to 26× coverage aiming for a genome of 300 megabases (Mb). The reads N50/N90 were estimated at 7,958/1,719.
Assembling reference genome
Reads aligning with the Mitochondrion genome were identified using BLASTN and filtered from the raw reads before assembling the genome. The remaining ONT reads were assembled using the Flye software (v.2.6), with –min-overlap 3,000 to increase stringency at the initial overlay step, and default parameters including five rounds of polishing through consensus, contigs were additionally polished two times with Medaka (v.0.11.2). Illumina paired-end RNA dataset is assembled using CLC Assembler (CLC Assembly Cell). Both RNA assemblies and paired-end 10X genomic dataset (unpublished data) were mapped onto the contigs using minimap2 (v.2.16) and BWA mapper (v.0.7.17), respectively, and the assembly was further polished using PILON (v.1.20) to error correct potential low-quality regions. The resulting assembly yielded a genome of 307 Mb, assembled into 1,533 contigs ranging from 10,840,357 to 100 basepairs (bp) and an assembly-N50 of 1.670 Mb. Moreover, the BUSCO completeness analysis using the Arachnida (odb10) reference set confirmed our assembly represents the complete genome C:94.8%(S:89.1%,D:5.7%),F:0.9%,M:4.3%,n:2934 (=arachnida_odb10), only missing 126 genes from the whole reference set. Knowing that BUSCO only gives a rough estimation, we remain confident that this assembly represents well the bulb mite genome.
Flow-cytometry
Whole individual R. robini were homogenized in 500 μl of ice-cold LB01 detergent buffer along with the head of a male Drosophila melanogaster (1 C = 0.18 pg) as an internal standard. The homogenized tissue was filtered through a 30-μm nylon filter. Then 12 μl of propidium iodide with 2 μl of RNase was added, and stained for 1 h on ice in the dark. All samples were run on an FC500 flow cytometer (Beckman-Coulter) using a 488-nm blue laser, providing output as single-parameter histograms showing relative fluorescence between the standard nuclei and the R. robini nuclei. Six replicate samples were run to account for variation in fluorescence outputs. The genome size of R. robini was estimated at 0.30 pg, or about 293 Mb, and consistent with estimation of size from the genome assembly described above.
Mitochondrial genome
ONT reads aligned with R. robini mitochondrion genome were de novo assembled with Flye (v.2.6) assembler and polished with Racon. Mitochondrion genome is assembled in one single contig with a size of 15,335 bases.
Gene prediction
On the polished final genome, protein coding genes have been predicted. For this, AUGUSTUS was used including hints coming from R. robini RNA-sequencing (RNA-seq) (samples SRR3934324, SRR3934325, SRR3934326, SRR3934327, SRR3934328, SRR3934329, SRR3934330, SRR3934331, SRR3934332, SRR3934333, SRR3934335, SRR3934337, SRR3934338 and SRR3934339 from the PRJNA330592 BioProject deposited at the National Center for Biotechnology Information (NCBI) Short Read Archive) and proteins coming from highly curated Tetranychus urticae (v.2020-03-20) as well as proteins from the previous version of the unpublished, Illumina-sequenced R. robini genome (https://public-docs.crg.eu/rguigo/Data/fcamara/bulbmite.v4a/). The PE RNA-seq reads were mapped on the genome using HISAT2 (-k 1 —no-unal) and further processed with Regtools to extract junction hints and filtered for junctions with a minimum coverage of 10. All the RNA-seq reads were also assembled with CLC Assembly Cell (v.5.2.0) software, setting the word size for the Bruijn graph at 50 and maximum bubble size at 31. The reads were assembled into 689,563 contigs (ranging from 10,675 to 180 bp), which were later mapped on the genome with GenomeTheader to generate complementary DNA hints. Protein hints were generated by using with Exonerate (v.2.2) with Protein2Genome model. To reduce the amount of overprediction due to repeated elements (transposable elements, simple sequence repeats) we de novo predicted high abundant repeats using RepeatModeler. The accompanying parameter file for extrinsic data for AUGUSTUS was adapted to include these hints as well as the softmasking of the genomic sequence. The resulting gene predictions from AUGUSTUS were further curated with EvidenceModeler using the same extrinsic data. The BUSCO analysis confirmed that our gene prediction indeed captured the expected genes well (C:94.6%(S:86.3%,D:8.3%),F:0.4%,M:5.0%,n:2934 (=arachnida_odb10)). The final predicted gene set was subsequently processed to be uploaded into ORCAE (https://bioinformatics.psb.ugent.be/orcae/overview/Rhrob)98.
Resequencing
Genomic sampling and mapping
For genomic analyses we sampled material from each of the morph selection lines (n = 8) at F1, F12 and F29. Following the experimental evolution protocols, after the first 24 h of egg laying all adults were transferred to a new container (described above) for a second 24 h to lay eggs and from these second dishes genomic material was sampled. On maturation, adults were transferred to and kept for 3 days in containers. Adults were then randomly selected and placed into Digestion Solution for MagJET gDNA Kit (F1&12) or ATL buffer (F29) before freezing at −20 °C. From each population two samples were collected consisting of 100 individuals (1 × 100 females and 1 × 100 males of random morph), the two samples separated by sex were used as technical replicates. The tissue from the 100 individuals within each sample was homogenized and DNA was extracted by Proteinase K digestion (24 h) followed by standard procedures using MagJET Genomic DNA Kit (ThermoScientific, F1&12) or DNeasy Blood and Tissue (Qiagen, F29). DNA concentration was controlled with the Qubit double-stranded DNA HS Assay Kit and DNA quality was assessed on agarose gels. The library preparation was performed using NEBNext Ultra II FS DNA Library Prep kit for Illumina.
Whole-genome resequencing was carried out by National Genomics Infrastructure (Uppsala, Sweden) using the Illumina Nova-Seq 6000 platform with S4 flow cell to produce 2 × 150 bp reads (average 160.7 × 106: range 130.7 × 106 − 189.9 × 106). Adaptors were trimmed from reads using Trimmomatic99 software (v.0.39) and unpaired reads discarded. Fastq files were mapped to the assembled genome with bwa mem100 (v.0.7.17-r1188) using default settings. Sam files were converted to bam files, sorted, duplicates marked and ambiguously mapped reads removed using samtools101 (v.1.9). On average, 90% (range, 86–93%) of the reads from each sample were mapped successfully, of which an average of 17% (range, 15–19%) were marked as duplicates. This left us with an average of 117.7 × 106 pair end reads per sample, ranging between 99.6 × 106 and 145.9 × 106 (Supplementary Table 1).
Genomic analysis
File preparation and filtering
Preparation of files used in genomic analysis was done as follows: bam files were converted to a pile-up file using samtools, following which indels and surrounding windows (5 bp either side) were filtered, using identify-genomic-indel-regions.pl and filter-pileup-by-gtf.pl in PoPoolation102 (v.1.2.2) to avoid false SNPs, with the resulting filtered pile-up file converted to a sync file using mpileup2sync.pl in PoPoolation2 (ref. 103) (v.1.201). Using custom python scripts, the distribution of coverage from each sample (single sex) was determined by recording the coverage of positions every 10 kb across the genome from the sync file to give information on expected coverage (Supplementary Fig. 1). On the basis of this, we filtered the sync and pile-up files to contain only regions within a range of informative coverage, where the mean coverage of all samples at every position was between 50% of the expected coverage and 200% of the expected coverage (56×, range 23−112×). The pile-up and sync files containing individual male and female samples (48 in total) were then merged by sex to give files containing allele frequencies from 24 samples (eight populations across three generations), each consisting of allele frequencies of 200 individuals (100 males and 100 females, above) and used in all subsequent analysis (unless stated otherwise). Similarly, we drew coverage of a position every 10 kb from each sample in the sex-merged sync file to determine a distribution from which we decided to subsample to (Supplementary Fig. 1). We putatively identified X-linked contigs (below) and excluded them autosomal analysis. A similar, but, separate analysis on genes and SNPs from X-linked contigs was performed by using different parameters (below).
Estimating nucleotide diversity
Using PoPoolation we determined various estimates of genetic diversity per sample (that is, 24 sex-merged samples). The pile-up file from each sample was subsampled using subsample-pileup.pl to a coverage of 63× (max coverage, 252×) to standardize estimations of genetic diversity across the genome, between populations and across generations. First, nucleotide diversity (Tajima’s Pi, π) and number of segregating sites (Watterson’s theta, ϴ) were estimated within genes. We performed analysis of exons using Syn-nonsyn-at-position.pl, in which genetic diversity of synonymous and non-synonymous positions were determined. Further analysis of overall genetic diversity within exons and introns were performed using Variance-at-position.pl, Tajima’s D (D) also estimated in the former. We used a minimum count of three (equal to a minor allele frequency of roughly 5%) for a SNP to be called, and a phred score >30 and a pool size of 400. Further analysis using 10 kb sliding windows (step size 10 kb) across the genome were performed using Variance-sliding.pl, and also included estimation of D. Estimates of D require the minimum count to be 2, but otherwise all the same parameters were used.
We filtered genes to be included in our analysis (and all subsequent analysis) on the basis of a number of criteria. On the basis of extensive RNA-seq data from both males and females (Plesnar-Bielak, unpublished data with NCBI accession number PRJNA796800), we only included genes in our analyses that were expressed at a mean level of fragments per kilobase of transcript per million mapped reads >1 across 72 samples originated from both sexes and both morphs rearing in three different temperatures (18, 23 and 28 °C). A further filtering step was performed to remove genes with inconsistent mapping between samples, only genes with >60% exons mapped to (calculated from positions used to calculate parameters in the Syn-nonsyn-at-position.pl π outputs), with 63−252× coverage, in all 24 samples were included in the analysis. The final dataset contained 13,389 autosomal genes and subsequently used to filter other datasets to retain this set of genes only (see Supplementary Table 8 for a list of genes). Similarly, windows were discarded from outputs if <60% had been mapped to with 63−252× coverage; when comparing the estimation of genetic diversity between treatments or generations, every window in all 24 samples had to meet these criteria. This means that all comparisons are conservative and based on the same genes or windows, and therefore unlikely to be biased by any differences in mapping between samples.
As our data included estimates of genetic diversity from the same population across multiple generations, we performed analysis by repeated measures analysis of variance (ANOVA), which takes into account the non-independence of samples. Comparisons were made on the mean values of each experimental evolution population.
Estimating X-linked diversity
We then repeated the above analysis on X-linked contigs (below) using the same parameters unless stated otherwise, and also performing the SNP based analysis below. Initially we ran the analysis using 75% of the target coverage and pool size used for autosomes that is a target coverage of 47× and a maximum coverage of 189×, with a pool size of 300. Following filtering steps, it was clear that two samples (PS17 and PS21) with relatively low coverage (Supplementary Fig. 1 and Supplementary Table 1) were having a large effect on filtering steps (that is, >60% of genes being mapped to in all 24 samples) and reducing the final X-linked dataset to contain fewer than 200 genes. We therefore opted to reduce the target coverage further to 40×, in an attempt to retain more genes. This slight reduction of target coverage increased the number of genes in the final dataset substantially to 587 genes. We therefore opted to use a minimum coverage of 40× in all analysis of X-linked SNPs, genes and windows.
Diverging SNPs
To determine divergent SNPs between F- and S-lines, we extracted the allele frequencies of all samples from the sex combined sync file. Samples from F29 were then used to filter the entire dataset to only contain SNPs on the basis of a number of criteria. First, positions within all samples were required to have a coverage >63× and <252×. Second, the frequency of minor alleles (calculated as coverage − major allele count) from all samples combined had to be >5% (that is, the average of all samples but not necessarily above >5% in all samples). Thus, our dataset contained only positions with the target coverage in all F29 samples and in which polymorphisms were unlikely to be a consequence of sequencing errors. After this filtering we were left with roughly 6 million SNPs used in further analysis. We performed a GLM, at each position by comparing the count of the major allele against counts of minor alleles at F29, to determine consistent allele frequency changes between treatments70. If any population had minor or major allele count of 0, +1 was added to minor and major alleles from all samples. To correct for multiple testing, we converted P values to q values using the qvalues R package (v.2.14.1)104 and applied a FDR with a q < 0.05, leaving 24,189 consistently diverged SNPs. Of those SNPs that we classified as diverged at F29, we then performed GLMs at these positions on F1 samples to examine whether they began the experiment diverged.
Initial SNP frequencies
We then compared the initial frequency of alleles (that is, at F1) of diverged SNPs to the genomic average (autosomes only). From the F1 samples, we determined the average allele frequency of all populations (that is, treating all F1 samples as a panmictic population), we then randomly sampled 24,189 positions and recorded the median frequency of minor alleles; this was repeated 10,000 times to draw a distribution. From this distribution, we then determined CIs to examine whether the median frequency of minor alleles of diverged SNPs differed from the genomic average.
Position of diverged SNPs
Using bedtools105 (v.2.27.1) we determined which genes (exons) contained diverged SNPs. Next also using bedtools, the genome was split into 10 kb windows and we determined windows that contained at least one significantly diverged SNP. By then drawing 24,189 random positions from all possible SNPs and counting the number of autosomal windows at least one SNP was within, and repeating this 10,000 times, we were able to draw a distribution and determine CIs of the number of windows containing random SNPs. This was used to examine whether diverged SNPs were distributed randomly across the genome. Position of diverged SNPs were visualized by Manhattan plots using the R package qqman106.
Regions or genes containing diverged SNPs
The ratio of non-synonymous (PN) to synonymous (PS) segregating sites (PN/PS) was compared between exons that contained diverged SNPs against those that did not by Wilcoxon signed-rank test on the basis of the average across replicates of F1 samples. Additionally, to test a specific hypothesis that genes containing diverged SNPs that were fixed in F-lines are under stronger purifying selection than genomic average, PN/PS of exons containing them were compared to those that did not contain diverged SNPs fixed in F-lines. Due to the relatively low numbers in the former set of genes, we compared these two sets using a random sampling approach, in which we randomly sampled 78 times from the set of genes that did not contain diverged SNPs, and the median was calculated. This was repeated 10,000 times to draw a distribution and calculate 95% CIs for the median.
Furthermore, genetic diversity was compared between genes (exons) and 10 kb windows that contained diverged SNPs against those that did not. For this purpose, we calculated the mean π of F1 samples for both exonic diversity and that within 10 kb windows, and compared the two groups using Wilcoxon signed-rank tests due to non-normal distributions.
Finally, we explored the top 10% genes and 10 kb windows with highest values of D (that is, those with signatures of balancing selection) calculated as an average across F1 samples. We investigated whether the top 10% genes and windows were enriched for those containing diverged SNPs using chi-square test analyses. We similarly tested whether the top 10% windows were enriched for those containing SNPs fixed in S-lines, to test the specific hypothesis that these regions were initially under balancing selection due to sexual antagonism and the SNPs were subsequently lost in those lines when sexual antagonism was relaxed. Finally, to test whether the top 10% D set at F1 was more likely to maintain high values of D across the experimental evolution in F-lines compared to S-lines, we compared mean D values of this set of 10 kb windows between treatments at F29 using a simple t-test.
Identifying X-linked contigs
As with most other Acarid mites107, R. robini has a XO sex determination system, with males being the heterogametic sex (Supplementary Fig. 4). On the basis of predicted differences in read coverages of contigs between male and female samples we identified putative X-linked contigs using the Illumina reads from all experimental evolution populations. We calculated the mean coverage of each contig for individual male and female samples, with contig coverage standardized by the mean overall sample coverage. Autosomal contigs are expected to have a ratio between female and male samples of 1:1, whereas X-linked contigs are expected to have a ratio of 1:0.5. These expected differences in the latter are due to females having two copies of the X chromosome but males only having a single copy. A similar approach was used in Callosobruchus maculatus57 to putatively assign contigs as either autosomal or as a sex chromosome (X- and Y-linked). Inspection of Supplementary Fig. 5 shows that few contigs conform to the expectation of X-linked contigs, but that a number of contigs do have relatively low male coverage. It appears that on average male samples have slightly higher than expected coverage; this is probably a consequence of males only having a single sex chromosome and therefore an excess of autosomes and X-chromosome reads are found in male libraries compared to female libraries. With this in mind, we drew an arbitrary cut-off in the ratio of coverage between female and male contigs of 1:0.75, with contigs with a ratio above this being assigned as autosomal, and those with a ratio below this assigned as X-linked contigs (Supplementary Fig. 5).
Estimating N
e, selection coefficients and simulating the impact of drift with different N
e
Drift under different N
e
An estimation of Ne for each population was determined by the changes in allele frequency of all autosomal SNPs between F1 and F29 using the R package poolSeq108 (v.0.3.5) and the estimateNe() function. As these results showed a difference in Ne between morph selection treatments, we performed simulations to determine whether this differences in Ne and consequently changes in drift could explain our patterns of diverged SNPs. Using poolSeq and the wf.traj() function we simulated two populations with Ne the same as our mean estimation of Ne of the F-lines (Ne = 370) and S-lines (Ne = 460), with generation of sampling (excluding F12), census size and sample size identical to those used in the experiment. To save computing time, we ran 20 identical simulations and combined output files. Each simulated population consisted of 1.5 million independent SNPs evolving in a Wright–Fisher population, the starting frequencies of minor alleles were randomly drawn between the range of 0.01 and 0.5. From these two simulated populations we sampled 50,000 times to determine whether allele frequencies consistently diverged. We approximated initial frequency of SNPs in the following way: by binning initial allele frequency from genetic F1 SNP data (<1, 1–5, 5–10% and so on in steps of 5%) we calculated a proportion of total SNPs within each bin. Using this proportion, we sampled from each simulated bin from both simulated populations in proportions matching those from genetic data (that is, 50,000 × proportion SNPs in bin). Each sample consisted of drawing four positions (without replacement) from each of the simulated populations from the same F1 bin, thus, adding a small amount of noise to the initial allele frequency of each sample. To add some further noise associated with differences in coverage, for each sample we drew a random number from a Poisson distribution resembling our actual target coverage (lambda, 125), by then drawing eight random numbers between −20 and +20 (approximate variation estimated from genetic data) and adding these to the random number from the Poisson distribution. Each of the simulated proportions of allele frequencies at F29 in the sample were then multiplied by one of these eight numbers at random, and rounded to the nearest integer, to gain counts of minor and major alleles and introduce some variation in coverage comparable to our molecular data. As with molecular data, we discarded simulated positions that did not meet our criteria of having an average proportion of minor alleles from all simulated populations >5% and if the coverage from any simulated position was <63× and >252×. From those samples that remained (>900,000), GLMs were performed (identical to above) on the simulated major and minor allele counts. Using a FDR with a q < 0.05 no simulated positions reached significance.
Selection coefficients
Next using the estimateSH() function in poolSeq we estimated selection coefficients (s) within each morph selection treatment on the set of diverged SNPs. Initially, we tried to allow dominance to be accounted for at each position using method = ’automatic’, but due to low number of time points at which sampling was performed this was not possible and all outputs reverted to using codominance. We therefore opted to use fixed codominance (h = 0.5), which also enabled us to calculate a P value by running 1,000 simulations for each position to estimate if s differed significantly from drift (that is, s = 0) within each morph selection treatment. Reliable estimates of s are not feasible if allele frequencies are too low, therefore, we were unable to determine s for every position in each morph selection treatments. To account for multiple testing, we applied a FDR with a q < 0.05.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
