Within-population sperm competition intensity does not predict asymmetry in conpopulation sperm precedence

Postcopulatory sexual selection can generate evolutionary arms races between the sexes resulting in the rapid coevolution of reproductive phenotypes. As traits affecting fertilization success diverge between populations, postmating prezygotic (PMPZ) barriers to gene flow may evolve. Conspecific sperm precedence is a form of PMPZ isolation thought to evolve early during speciation yet has mostly been studied between species. Here, we show conpopulation sperm precedence (CpSP) between Drosophila montana populations. Using Pool-seq genomic data we estimate divergence times and ask whether PMPZ isolation evolved in the face of gene flow. We find models incorporating gene flow fit the data best indicating populations experienced considerable gene flow during divergence. We find CpSP is asymmetric and mirrors asymmetry in non-competitive PMPZ isolation, suggesting these phenomena have a shared mechanism. However, we show asymmetry is unrelated to the strength of postcopulatory sexual selection acting within populations. We tested whether overlapping foreign and coevolved ejaculates within the female reproductive tract altered fertilization success but found no effect. Our results show that neither time since divergence nor sperm competitiveness predicts the strength of PMPZ isolation. We suggest that instead cryptic female choice or mutation-order divergence may drive divergence of postcopulatory phenotypes resulting in PMPZ isolation. This article is part of the theme issue ‘Fifty years of sperm competition’.

: Summary of prezygotic isolation between populations of Drosophila montana from Crested Butte, Colorado, USA, Oulanka, Finland, and Vancouver, Canada. Edge weights reflect absolute contributions of premating sexual isolation (dashed) and non-competitive PMPZ isolation (solid) to total reproductive isolation between each pair of populations (from Table 3 in Jennings et al. (2014)).

Population genetics DNA sequencing
Genomic DNA was extracted from pooled samples of 50 females per population and sequenced using an Illumina HiSeq 2000 to produce paired-end reads (90 + 90 bp, insert size = 500 bp) at the Beijing Genomics Institute (Parker et al. 2018). Reads were trimmed, filtered and down-sampled to ensure equal representation of each population (Parker et al. 2018). Sequencing was performed in 2013, 5 years after collection for Oulanka and Vancouver, and 6 years for Colorado (around 50 generations).

Demographic modelling
All single nucleotide polymorphisms (SNPs) were filtered using a minimum coverage of 30 and we randomly subsampled SNPs to end with a uniform coverage of 30. To ensure SNPs included in the analysis were unlikely to be linked, we then randomly down-sampled based on how far SNPs were apart by selecting 1 in every 100 SNPs.
We indexed the D. montana Illumina reference genome and mapped Pool-seq reads using BWA-MEM v.0.75 (Li 2013). We sorted reads, removed duplicates, and converted mapped reads using SAMtools v.1.6 (Li et al. 2009) and Popoolation2 (Kofler et al. 2011). All single nucleotide polymorphisms (SNPs) were filtered using a minimum coverage of 30 and we randomly subsampled SNPs to end with a uniform coverage of 30 before converting read counts into a format suitable for δaδi (Diffusion Approximation for Demographic Inference. To satisfy the assumption of independence between SNPs, we down-sampled from 860,228 to 8,000 SNPs by selecting 1 in every 100 SNPs post-filtering to ensure SNPs in close proximity were not used in the demographic analysis. We tested seven demographic models for the history of divergence between D. montana populations from Crested Butte, Colorado, USA, Oulanka, Finland, and Vancouver, Canada (Fig. S2).
(1) a no-migration model where migration was not included; (2) a full-symmetrical migration model where migration was symmetrical and present between all populations; (3) an adjacent-migration model where migration was possible between Oulanka and Vancouver or between Colorado and Vancouver, but not between Colorado and Oulanka. Additionally, we included models incorporating periods of isolation falling under two categories: ancient migration models and secondary contact models. Ancient migration models included (4) a model assuming migration ceased at the time of the Colorado-Vancouver split, and (5) a model assuming migration persisted after the Finnish-North American split and the Colorado-Vancouver split but ceased shortly after. For secondary contact models, we included (6) a model with isolation between populations with a recent secondary contact event between Colorado and Vancouver after divergence, and (7) a model assuming migration during the Colorado-Vancouver divergence, with a period of secondary contact beginning after the appearance of both Colorado and Vancouver populations. Migration is assumed to be symmetric for all models.
We fitted all seven models to the 3D-AFS and performed 3 rounds of model optimizations using the Nelder-Mead method (Nelder and Mead 1965). In the first round, minimum (0.01) and maximum (10) parameters were 3-fold perturbed with 10 replicates and a maximum of 5 iterations. The best fitting parameters from the first round were used as starting parameters for the second round, in which parameters were 2-fold perturbed with 10 replicates and a maximum of 5 iterations. Finally, the best fit replicate from the second round was used to produce starting parameters for the third round, in which we ran 15 replicates with parameters perturbed 1-fold with a maximum of 5 iterations.
We estimated divergence times by first solving the equation for θ, the population mutation rate, to derive N ref , the reference effective population size: where µ is the mutation rate for D. melanogaster (2.8x10 −9 ) (Keightley et al. 2014) and L, the effective length was calculated by multiplying the total number of SNPs aligned to the reference genome and converted to sync format (31,800,188), by the number of SNPs that entered the analysis divided by the number of SNPs remaining before down-sampling (8,000/860,228). We then converted parameter estimates to real biological units of time, τ i, using the equation: where G, the generation time, was set as one per year based on observations from Aspi et al. (1993) and Throckmorton (1982).
Using equations (1) and (2) we estimated divergence from the ancestral population to have occurred 1,743,164 years ago (T1 + T2), and between Colorado and Vancouver shortly after (T2 = 1,710,355 years ago). We advise caution when interpreting these estimates however. Despite studies showing accurate allele-frequency estimation using Pool-seq (Futschik & Schlötterer 2010;Gautier et al. 2013), it is difficult to appreciate how the use of read count data derived from a DNA pool might affect demographic inference and divergence time estimates. Here, for example, we observe some evidence of an excess of rare variants in our data, though it is unclear whether this represents demographic and/or selective processes, or inaccurate allele frequency estimation caused by inherent statistical challenges associated with Pool-seq. To account for this, we stringently filtered and down-sampled our dataset to retain only reliable sites. Nonetheless, an excess of rare variants may have biased the demographic inference performed here. Additionally, the generation time (G), effect seqeunce length (L) and mutation rate (µ) used, whilst these are reasonable assumptions, may have also contributed to over-or underestimation of divergence times.  Figure S2: Graphical representation of demographic models tested. mA, ancient migration between Oulanka and precursor North American population; m1, ancient migration between Vancouver and Oulanka; m2, ancient migration between Vancouver and Colorado; m3, ancient migration between Oulanka and Colorado; Initial and secondary divergence are represented by T1(a) and T2, where T1(b) and T3 correspond to the start or end of periods of isolation. Population splits are denoted by T1 (between Oulanka and North America) and T2 (between Vancouver and Colorado). T1(a) and T1(b) denote the start or end of periods of migration following the split between Oulanka and North American populations. T3 denotes start or end of periods of migration following the split between Vancouver and Colorado.

Estimating within population genetic diversity
We mapped reads for each population to an updated D. montana PacBio reference genome supplied by colleagues (Poikela et al. in prep.) as above using BWA-MEM and conversion to pileup using SAMtools. To estimate θ W atterson , θ π and Tajima's D, we used NPStat v.1 (Ferretti et al. 2013) and calculated summary statistics in 1kb windows, with minimum coverage (10), maximum coverage (100), minimum base quality (30) and minimum allele count (2) filters. We averaged windows for each of 53 scaffolds (representing 95.5% of the genome) to obtain scaffold-wide estimates of summary statistics and evaluated differences between populations using Kruskal-Wallis rank sum tests followed by pairwise Wilcox rank sum tests corrected for multiple testing using the Benjamini-Hochberg method.  Figure S4: Genome-wide estimates of Tajima's D, pi (π), and Watterson's theta (θ W ) calculated in 1kb windows. Points show the mean for each scaffold, n = 53 in each case. Pairwise differences between populations were evaluated using Wilcox rank sum tests corrected for multiple testing using the Benjamini-Hochberg method: * p < 0.05, ** p < 0.01, *** p < 0.001.

Measures of postmating prezygotic reproductive isolation
We measured PMPZ isolation between Colorado and Vancouver using hatching success rates (a proxy for fertilisation success) (Jennings et al. 2014;Garlovsky and Snook 2018) in single and double matings, and assessed CpSP by estimating male offensive paternity share (P2) using the irradiated male technique (Boorman and Parker 1976). We performed experiments in 3 blocks for each focal female population separately, in which virgin females (n = 10 per cross-type per treatment per block) were presented with a virgin male from Colorado or Vancouver and allowed a 4-hour mating opportunity. After mating, males were discarded and females were transferred to a chamber in an oviposition manifold, which allows easy counting of numbers of eggs oviposited and assessment of hatching success (Jennings et al. 2014;Garlovsky and Snook 2018).
Manifolds were returned to the incubator for oviposition. Two days later, females were removed from their chambers, and placed in a vial with a second virgin male from Colorado or Vancouver and allowed a 4-hour remating opportunity. Remating females were returned to their chambers, placed back into the incubator and allowed to oviposit for a further two days on a new oviposition plate. We ensured no females mated more than once during either the first or second observation period. To measure hatching success the numbers of eggs laid on each oviposition plate were counted immediately after removing the female (for remating or discarding). The oviposition plate was then returned to the incubator for a further two days after which the numbers of unhatched eggs was counted again.
We assessed differences between cross-types in hatching success rates after the first and second mating (for females mated to fertile males only) using generalized linear models (GLMs) with binomial errors and a logit link, where the numbers of hatched ("successes") and unhatched ("failures") eggs was the response variable and the only predictor was the cross-type. We assessed differences in P2 using GLMs with binomial errors and a logit link where the binary response was the numbers of offspring sired by the second male ("successes") and numbers of offspring sired by the first male ("failures") (see supplementary material). Models included cross-type, irradiation order and their interaction as predictors. Females that did not remate, died prior to remating, or laid less than 10 eggs were excluded. Due to the low fertility in crosses between Colorado females and Vancouver males, we could not estimate P2 in the CVV cross, which was subsequently excluded from analyses (see supplementary material). We analysed responses in Colorado and Vancouver females separately as populations were never tested together. All statistical analyses were performed in R v.3.5.1 (R Core Team 2018). We used quasi-errors for GLMs where preliminary model inspection indicated overdispersion. Where appropriate we performed post-hoc Tukey's honest significant difference (HSD) tests using the glht function from the 'multcomp' package (Hothorn et al. 2008).  Figure S6: Hatching success (% eggs laid that hatched) after the second mating for females that mated with nonirradiated males. Points are observations and size show the number of eggs laid (i.e. weights). Cross-type denotes female population followed by first and second male. C, Colorado; V, Vancouver.

Con-population sperm precedence (CpSP)
We assessed male offensive paternity share (P 2 ) using the irradiated male technique (Boorman and Parker 1976). Virgin males were irradiated with a 100Gy dose of gamma radiation (dose rate 189.2 rads min −1 , 137 Cs source), rendering males 100% sterile, less than 24 hours before mating and housed individually overnight in food vials. To assess the level of fertility after a double mating, 'control' treatments consisted of mating females with two nonirradiated males in all possible crossing combinations with males from Colorado and/or Vancouver. We also included additional data for these control crosses collected in a pilot study (n = 30 per cross-type). To assess the efficacy of the irradiation technique females were mated with two irradiated males, both from Colorado or Vancouver. Experimental treatments consisted of all possible crossing combinations between females and males from Colorado and Vancouver, with either the first or second male irradiated. The proportion of eggs fertilised by the irradiated male after the second mating, P R , was calculated using equation (1) from Boorman and Parker (1976): where x is the observed proportion of developing eggs after a double mating from the second oviposition plate, p is the level of fertility observed in a double mating for a given 'control' cross-type, and z is the level of fertility observed in the double irradiated cross-type. For our calculation of P R we chose the maximum value of p observed (rather than the mean) to capture the full potential fertility of between-population crosses. As z = 0 in this case, the equation can be simplified to P R = 1 − x/p. Therefore, if the irradiated male mates first, then the proportion of eggs fertilised by the second male, P 2 = x/p. If the irradiated male mates second then P 2 = P R (Boorman and Parker 1976). The total number of eggs laid by each female after the second mating was multiplied by the calculated P 2 value and rounded to a whole number to give the estimated number of offspring sired by the second male. The remaining number of eggs expected to hatch were assigned to the first male.
As few eggs laid hatch when Colorado females mate with Vancouver males (CV hatching success = 0.084 ± 0.01 [mean ± standard error], n = 79, Fig. S5; CVV hatching success = 0.115 ± 0.033, n = 30, Fig. S6), we could not estimate differences in P 2 for the CVV cross-type which was subsequently excluded from P 2 analyses (Fig. S7). In Colorado there was a significant effect of irradiation treatment (F 1,82 = 20.41, p < 0.001), and the cross-type x irradiation interaction on P 2 (F 2,80 = 14.12, p < 0.001). Irradiated males mating in the second position had a greater P 2 than irradiated males mating in the first position in the CCV and CVC cross-types, but not the CCC cross-type (Fig. S7). In Vancouver the irradiation main effect was not significant (F 1,116 = 1.24, p = 0.267) but there was a significant cross-type x irradiation interaction (F 3,113 = 7.37, p < 0.001). In the VVV cross-type irradiated males in the second position had lower P 2 , whereas in the other Vancouver female cross-types irradiated males in the second position had greater P 2 (Fig. S7).

Interaction between coevolved and foreign male ejaculates in the female reproductive tract
We tested whether observed hatching success rates after a double mating in the nonirradiated 'control' crosses differed from the expected additive effect of two single matings, H total , using the following equation: where P 2 is the mean proportion of offspring sired by the second male in a given cross-type and H 1 and H 2 are the mean hatching success rates after a single mating for a female mated with a male from the first, and second, population denoted in a cross-type, respectively (Table 2). Hatching success (%) Figure S8: Effects of overlapping foreign and coevolved ejaculates on fertilisation success. White points show H obs , mean hatching succcess (% eggs laid that hatched) for females that mated with two nonirradiated males (data same as Fig. S6), red diamonds show H total , the expected additive hatching success of two single matings for a female mated with a male from the first and second population denoted in a cross-type, calculated given the estimated paternity share of each male (Table 2). Cross-type denotes female population followed by first and second male mate. C, Colorado; V, Vancouver.

Male relative reproductive investment
We analysed differences between populations in male relative reproductive mass as a measure of total ejaculate investment (Tomkins and Simmons 2002). Approximately 100 flies of mixed sex were collected from population cages and allowed to mate and oviposit for two days in plastic bottles covered with a molasses-agar oviposition plate and a drop of dried yeast paste. Controlled density vials (CDVs) containing food were seeded with 50 first instar larvae from oviposition plates. Adults were collected from CDVs on the day of eclosion and kept in single sex vials until sexually mature and then frozen at -20ºC. Population identity was replaced with a unique identifier prior to measurement. Males (n = 60 per population) were thawed and the entire reproductive tract (testes, accessory glands, ejaculatory duct and bulb) dissected on a pre-weighed piece of foil in a drop of dH20 which was then transferred to a second piece of pre-weighed foil. Carcasses and reproductive tracts were dried overnight at 60ºC before weighing (METTLER TOLEDO® UMX2 ultra-microbalance).
We analysed differences between populations in male relative reproductive investment using analysis of covariance (ANCOVA) (Tomkins and Simmons 2002). Log transformed dry reproductive tract mass was regressed against log transformed dry soma mass (body mass -reproductive tract mass) and the soma mass x population interaction. Colorado males weighed significantly less than Vancouver males on average (Welch Two Sample t-test, t = -9.84, df = 96.45, p < 0.001; Colorado = 0.58 ± 0.01 mg [mean ± standard error], n = 60, Vancouver = 0.77 ± 0.02 mg, n = 60; Fig. S9a). Colorado males reproductive tracts weighed significantly less than Vancouver males on average (Welch Two Sample t-test, t = -6.47, df = 114.14, p < 0.001; Colorado = 0.09 ± 0.002 mg, n = 60, Vancouver = 0.11 ± 0.002 mg, n = 60; Fig. S9c). The results of the ANCOVA revealed the population x log soma mass interaction was not significant (ANCOVA: F 1,116 = 0.995, p = 0.321) and so was dropped from the model. The reduced model showed a significant effect of population (ANCOVA: F 1,117 = 7.74, p = 0.006) and log soma mass (ANCOVA: F 1,117 = 8.63, p = 0.004) on reproductive mass. Log reproductive tract mass increased with log soma mass and Vancouver males had a higher intercept (Fig.  S9b). However, visual inspection of model diagnostic plots revealed four consistent outliers (Fig. A1). One Vancouver male had an unusually small body mass which resulted in high leverage (Fig. S9a). Three males (2 Colorado, 1 Vancouver) had relatively small reproductive tract masses (Fig. S9c). After removing the four outliers, Colorado males still weighed significantly less than Vancouver males on average (Welch Two Sample t-test, t = -11.36, df = 98.8, p < 0.001) and Colorado reproductive tracts weighed significantly less than Vancouver males on average (Welch Two Sample t-test, t = -6.91, df = 106.27, p < 0.001). Results of the updated ANCOVA model again revealed that the population x log soma mass interaction was not significant (ANCOVA: F 1,112 = 0.02, p = 0.884), and was subsequently dropped from the model. Visual inspection of model diagnostic plots showed no significant outliers (Fig. A2).

Male mating capacity
Virgin males (n = 20 per population) were mouth aspirated without anaesthesia into a vial containing two virgin females from their own population. Once copulation began unmated females were removed from the vial. The male was transferred to a new vial housing another two virgin females after mating ad libitum and we recorded the total number of matings that males performed within a 4-hour period. Mated females were returned to the incubator in their oviposition vial for 4 days, before being transferred to a second food vial for a further 4 days. The total number of adult flies emerging was counted for each female (combined across both oviposition vials).