Discordant Pleistocene population size histories in a guild of hymenopteran parasitoids

Signatures of past changes in population size have been detected in genome‐wide variation in many species. However, the causes of such demographic changes and the extent to which they are shared across co‐distributed species remain poorly understood. During Pleistocene glacial maxima, many temperate European species were confined to southern refugia. While vicariance and range expansion processes associated with glacial cycles have been widely documented, it is unclear whether refugial populations of co‐distributed species have experienced shared histories of population size change. We analyse whole‐genome sequence data to reconstruct and compare demographic histories during the Quaternary for Iberian refuge populations in a single ecological guild (seven species of chalcid parasitoid wasps associated with oak cynipid galls). For four of these species, we find support for large changes in effective population size (Ne) through the Pleistocene that coincide with major climate events. However, there is little evidence that the timing, direction and magnitude of demographic change are shared across species, suggesting that demographic histories in this guild are largely idiosyncratic, even at the scale of a single glacial refugium.


| INTRODUC TI ON
Natural populations change in size and range in response to biotic and abiotic factors over both ecological and evolutionary timescales. Given that sequence variation is the only source of information about the evolutionary past of many taxa, there is considerable interest in using genetic data to reconstruct the demographic history of populations (see Beichman et al., 2018, for a recent review).
A key question is to what extent past demographic changes have been shaped by climatic events and whether such responses are concordant across co-distributed taxa, or vary with species traits (Hickerson et al., 2010;Papadopoulou & Knowles, 2016;Stone et al., 2012).
A major aim of phylogeographic studies has been to understand the origins and determinants of spatial structure in natural populations (Avise et al., 1987;Hickerson et al., 2010). Systematic comparisons of demographic histories across co-distributed taxa offer the chance to identify shared abiotic triggers of demographic events (Bermingham & Avise, 1986;Hewitt, 2000;Hickerson et al., 2010). If divergence and admixture between populations are largely the outcome of shared climatic and geological history, one would expect co-distributed and ecologically similar species to show concordant histories. Concordance in vicariance times has indeed been found by a number of comparative phylogeographic studies (e.g. Leachè et al., 2007;Stone et al., 2012;Wan et al., 2021;Xue & Hickerson, 2020).
However, most species are likely to have undergone complex range shifts and, given that range changes are likely to depend on multiple traits, the idea of strict temporal concordance across species seems implausibly simplistic biologically (Papadopoulou & Knowles, 2016).
For example, many temperate plant and animal taxa (including humans) have colonized the Western Palaearctic in a series of westwards range expansions. Overlaid onto this older longitudinal colonization history is a cyclical history of latitudinal range changes (Hewitt, 2000;Hofreiter & Stewart, 2009). Starting with the onset of major glacial cycles in the mid Pleistocene, the ranges of temperate taxa in Europe were restricted to southern refugia during glacial maxima, primarily in Iberia, Italy and the Balkans (Feliner, 2011;Hewitt, 2000;Hofreiter & Stewart, 2009) and expanded northwards during interglacial periods. As a consequence, modern refugial populations commonly harbour higher genetic diversity than northern populations founded by postglacial range expansion (e.g. Hewitt, 2000;Rokas et al., 2003;Tison et al., 2014;Vitales et al., 2016) (but see Comps et al., 2001). While the structuring of genetic diversity by refugia is ubiquitous (Hewitt, 1996), recent comparative studies of the population and speciation history of European taxa have, perhaps unsurprisingly, shown no support for temporal concordance (Bunnefeld et al., 2018;Ebdon et al., 2021). Bunnefeld et al. (2018) used whole-genome sequence data to compare the timing of divergence and admixture events between three southern refugial populations across 13 co-distributed species in a tritrophic Western Palaearctic community comprising oak (Quercus) host plants, cynipid gall wasp herbivores and chalcid parasitoid natural enemies (Stone et al., 2002). The latter are obligate specialists of cynipid galls, allowing the guild of parasitoids to be considered in ecological isolation. Many of the component species show broad longitudinal distributions, extending from Iberia in the West to Iran in the East, and all species studied to date are genetically structured into major southern refuge populations (Lohse et al.,2010Nicholls et al.,2010Nicholls et al., , 2012Petit et al., 2002;Stone et al., 2012Stone et al., , 2017. Bunnefeld et al. (2018) showed that seven out of the 13 community members studied colonized refugia in the Balkans and Iberia through westwards range expansion from an eastern origin one or more glacial cycles in the past. However, four species had a discordant, Western origin and two showed no clear signal of a longitudinal history. Bunnefeld et al. (2018) also found little evidence for temporal concordance: although divergence and admixture between refugial populations are broadly concentrated in the late Pleistocene, the timing of divergence and admixture between refugia is to a large extent idiosyncratic across species. However, in order to explore a plausible space of necessarily complex demographic scenarios which capture aspects of both the longitudinal colonization history (divergence between refugia) and the more recent history of latitudinal range shifts (admixture between refugia), Bunnefeld et al. (2018) had to make drastic simplifying assumptions about the second major axis of demographic history: population size change. Specifically, their analyses assumed a single fixed effective population size per refugium. This immediately raises the question of whether the lack of concordance of species' histories found by Bunnefeld et al. (2018) is unique to the broad spatial and deep temporal scale captured by their models, or a more general property of the history of ecological communities which holds over a range of spatio-temporal scales. There are several reasons to expect the degree of temporal concordance to depend on spatial and temporal scales (Papadopoulou & Knowles, 2016). In particular, colonization and admixture histories may be more temporally discordant than population size changes, simply because the former depend on longdistance dispersal events which are both a function of a species' dispersal ability and inherently random (Smith et al., 2014).
Here, we use whole-genome sequence (WGS) data to quantify and compare the histories of population size change within a single glacial refugium across seven chalcid parasitoid species in the oak gall wasp community previously studied by Bunnefeld et al. (2018).
We focus on the Iberian refugium, which for most species in the oak gall wasp parasitoid community-including four out of the seven parasitoids studied here (Table 1)-represents the end point of a longitudinal expansion history: Iberian populations have lower genetic diversity than more eastern refugia and show little evidence of population structure (Nicholls et al., 2010;Rokas et al., 2003;Stone et al., TA B L E 1 Maximum-composite likelihood estimates (MCLE) of parameters under the single-step-change model for Iberian populations of seven species of chalcid parasitoid wasps under a model of a single-step change. Estimates for the time of N e change (T) are given in thousands of years ago (kya). 95% CIs are shown in brackets. The maximum N e is scaled relative to as Max [N e ]∕N e, . The second column gives the origin of the longitudinal expansion history as inferred by Bunnefeld et al. (2018)  Ormyrus pomaceus West 0.00230 207 (207-208) n/a n/a n/a Eurytoma brunniventris East 0.00569 523 (521-525) n/a n/a n/a Cecidostiba fungosa East 0.00258 234 (233-235) n/a n/a n/a 2007) (for an exception, see Rokas et al. (2001) from the Balkans in each species), whose haploid genome facilitates analysis (Bunnefeld et al., 2018;Hearn et al., 2014). WGS data allow for powerful demographic inference even for non-model organisms.
However, given fragmented reference genomes and limited sample sizes available for these species, inference methods must be chosen carefully and, ideally, should use linkage as well as allele frequency information. In particular, inference based only on the site frequency spectrum (SFS) (e.g. Gutenkunst et al., 2009) requires larger samples for accurate inference of recent population history.
We investigate the history of population size change in Iberia using two approaches that use the signal contained in genome-wide variation of a small sample of individuals in different ways (Bunnefeld et al., 2018;Li & Durbin, 2011;Lohse et al., 2011). We use a parametric maximum-composite likelihood (MCL) method based on a blockwise summary of sequence variation (Lohse et al., 2011) to fit a model of a single instantaneous step change in N e (hereafter referred to as the 'single-step-change analyses'). Using analytic likelihood calculations to fit such fully specified but minimally complex histories utilizes both allele frequency and linkage information and facilitates statistical comparison between species for the timing and magnitude of N e change. We also applied a non-parametric method commonly used to visualize population size change, the pairwise sequentially Markovian coalescent (PSMC) (Li & Durbin, 2011), which is based on the distribution of pairwise differences in minimal samples of two haploid genomes. PSMC generates a more resolved picture of N e change through time, including a qualitative assessment of the time at which populations have diverged. As the two methods make opposing assumptions, we will not attempt any formal comparison between them but rather compare inferences qualitatively.  Table S1).
The M. dorsalis individuals correspond to cryptic species 1 for this complex, as defined by Nicholls et al. (2010). For each species, we sampled five individuals from Iberia (the focal refugium) and one from Hungary (the Balkan refugium). Data for the Hungarian and two Iberian samples of each species were generated previously by Bunnefeld et al. (2018). We generated analogous data for the remaining three Iberian individuals using the protocols described by

| Inferring single-step-change models in population size
We fitted a model of a single, instantaneous step change in N e using the framework for blockwise likelihood calculations developed by Lohse et al. (2011). The single-step-change model includes three parameters: the scaled mutation rate = 4N 0 l (where N 0 is the current N e and l is block length); T, the time of N e change measured in 2N e generations; and = N 0 ∕N 1 , the relative magnitude of the N e change (where N 1 denotes the N e prior to the step change).
Following Bunnefeld et al. (2015), we summarized sequence variation in short blocks of a fixed length l by the (folded) blockwise site frequency spectrum (bSFS) using the script block_cutter_vcf.py tion NMaximise (Data S1).
To generate blockwise data for each species, we applied the same quality filters used for calling SNPs to all sites; that is, we identified regions of the genome with a base quality >10 and mapping quality >20 in each individual from bam files by the CallableLoci walker of GATK (v3.4) (Van der Auwera et al., 2013). Only regions meeting these criteria in all five Iberian individuals were included in further analyses. Custom scripts were used to partition the data into blocks of a fixed length l of callable sites. l was chosen to be inversely proportional to the pairwise genetic diversity of each species (Table 2) such that blocks, each capturing a random configuration of closely linked variants had an average two pairwise differences between samples. This ensures that the information content per block, as well as any effect of intra-block recombination, is consistent across species. Blocks with a physical span (including non-callable sites) of >2l or more than five 'None' (uncalled) sites were removed.
We assessed support for a step change in N e relative to the  (Table S3). Each data set had the same total length as the real data (after filtering) and was partitioned into 5,000 windows of sequences (for computational efficiency). Both a null model of constant N e and a step-change model were fitted to each replicate, and the 95% quantile of the difference in support between models was compared to that of the real data. Confidence intervals (CIs) of parameter estimates were obtained via an analogous parametric bootstrapping procedure: we simulated 100 data sets with recombination under the inferred step-change model and fitted that model to each simulation replicate. CIs were obtained as 2.5% and 97.5% quantiles of the distribution of parameter estimates and were centred around the point estimates of parameters obtained from the real data.

| Reconstructing ancestral population size with PSMC
The pairwise sequentially Markovian coalescent (PSMC) (Li & Durbin, 2011) was used to infer a history of population size change for each species. PSMC is a non-parametric method that reconstructs a tra- smaller block lengths ( Figure S1). We chose a block length of 25 bp for all analyses which minimizes these biases without excessive computational cost.
We inferred 30 free interval parameters across 64 time intervals (with the option -p '28*2+3+5'). The maximum coalescence time (−t) was set to 5, the initial value of ∕ (-r) to 1. 100 bootstraps were performed for each run.  (Table 2) as N e, = ∕(4 ) and computed the following measure of N e change: Max[N e ]∕N e, . Unlike the size of the step change ( ), this measures the maximum N e relative to the average over the entire history of a sample and is therefore expected to be smaller than .

| Calibrating the timing of events
Time estimates were converted into years assuming a mutation rate of 3.46 × 10 − 9 mutations per base per generation estimated for Drosophila melanogaster (Keightley et al., 2009). We assumed two generations per year for all species with the exception of M. stigmatizans which has a single generation per year . While this calibration allows comparisons with the estimates obtained by Bunnefeld et al. (2018) (calibrated in the same way), these absolute times are likely underestimates and should be treated with caution (given that we use a spontaneous mutation rate estimate and assume selective neutrality). We stress, however, that comparing the relative timing of demographic events across this set of parasitoids only relies on the assumption of the same per generation mutation rate across species rather than any absolute calibration.

| Cross-population PSMC analyses
To test whether population size changes in each species occurred before or after the divergence between the Iberian and Balkan refuge populations, we compared the PSMC trajectories of Iberian pairs and cross-population pairs (one Iberian and one Hungarian individual).
Any divergence between Iberian and the Balkan refuge populations should be visible as a separation between within-population (Iberian pair) and cross-population (Iberia-Balkan pair) PSMC trajectories.
Specifically, in the absence of significant post-divergence gene flow, we expect N e estimates for the cross-population pairs to increase from the time of the population split due to the accumulation of genetic differences between the populations. Thus, population size changes that occurred in an ancestral population (potentially outside Iberia) should predate the divergence of within and cross-population trajectories. In contrast, we would expect demographic events unique to the Iberian population to happen after divergence of the within-and cross-population PSMC trajectories. We also compared our PSMC divergence time estimates to those made by Bunnefeld et al. (2018) under explicit models of divergence and admixture.

| Potential population structure in Iberia
Population structure within any assumed panmictic population can confound inferences of past demography (Bunnefeld et al., 2015;Gattepaille et al., 2013). In particular, N e estimates may be inflated due to divergence between demes (Gattepaille et al., 2013).
Likewise, when samples are taken from the same deme, the (structured) coalescent generates a mixture of very recent within-deme ancestry and older ancestry resulting from migration between demes (Wakeley, 2008), which can mimic the signatures of a bottleneck Bunnefeld et al. (2015). To test for population structure within Iberia, we repeated PSMC analyses on every pairwise combination of our five individuals. In a structured population, and given that the sampling locations were widely spaced across Iberia (Table S1), N e trajectories are expected to differ between pairs of haploid males sampled from the same versus different sub-populations (analogous to the within-and between-population analyses involving Iberian and Hungarian samples).

| Sensitivity analyses
The assumption of no recombination within blocks (required to make the composite likelihood estimation of the single-step-change model tractable) potentially biases parameter estimates. Specifically, recombination may lead to a downward bias of N e estimates (Bunnefeld et al., 2015;Wall, 2003). We used simulations in msprime (Kelleher et al., 2016) to quantify the effect of recombination on parameter estimates (Table S2) It is well known that both selective sweeps (Smith & Haigh, 1974) and background selection (Charlesworth et al., 1993) affect variation at neutral sites in the genome which, in turn, can bias demographic inference (Ewing & Jensen, 2016;Schrider et al., 2016). To explore the effect of selection on the single-step-change analyses, we fitted step-change models separately to blockwise data generated from genic and intergenic regions of the O. nitidulus genome. Genes were predicted ab initio with Augustus (Stanke & Morgenstern, 2005).
Nasonia vitripennis, a chalcid parasitoid (family Pteromalidae), was used as a training set. If selection has a strong effect on demographic inference, estimates of are expected to be lower for genic compared to intergenic regions, as most selection tends to reduce diversity both at selective targets and linked regions of the genome (Smith & Haigh, 1974). Similarly, we would expect estimates for the time of the step change, T, to be biased towards the present in genic compared to intergenic regions.

| RE SULTS
After filtering, the blockwise data sets ranged in total length from 54 to 151 Mb (Table S3). The pairwise alignments used as input for PSMC spanned a total length of 161 to 379 Mb (Table S3). Despite this difference in overall length, which is mainly due to missing data and the difference in sample size (two versus five individuals), estimates of average pairwise diversity, as measured by (Nei, 1972), agreed well between the two data sets (Table 1). Across species, estimates spanned an order of magnitude, from 0.00067 in M.

| Large changes in N e detected in four species
Taking nitidulus) and an increase in one species (T. auratus) ( > 1) ( Figure 1 and Table 1). In all four species, the change in N e in the PSMC trajectory ( Figure 1)  in N e is supported in neither species, and PSMC suggests comparatively small population size changes for both (Table 1 and Table 2).
However, while we have resisted the temptation to attempt any formal comparison between PSMC trajectories and single-stepchange analyses which would necessary be post hoc and potentially misleading (given the opposing assumptions of these methods, see  Table 1 and Table 2). Second, in one species, E. brunniventris, singlestep-change and PSMC analyses are hard to reconcile: we found no significant support for a step change in this species, that is, we cannot reject a null model of a single fixed N e . Yet, the PMSC trajectory indicates a substantial (but gradual) increase in N e (Table 2).
Additional analyses (see Discussion) suggest that demographic inferences for E. brunniventris may be affected by genetic structure and/ or the low contiguity of its reference genome.

| No support for temporal concordance of population size change
The four species with support for a step change in N e (M.  (Figure 2 and Table 2).

| Changes in N e occur after the divergence of Iberian populations
In five species, the PSMC trajectories of within-population (Iberia) pairs diverge clearly from the cross-population (Iberia versus Hungarian) trajectories (Figure 3). In contrast, little divergence between within and cross-population PSMC trajectories is visible for E. brunniventris and none for T. auratus (Figure 3). In all three species for which the single-step-change analyses support a decline in population size, the inferred time of N e change is more recent than the divergence of within-and cross-population PSMC trajectories ( Figure 3 and Table 2) and so must have occurred after Iberian and In both cases, cross-population PSMC trajectories decrease after divergence between the Iberian and Hungarian populations, suggesting that these populations have been connected by a period of gene flow immediately following divergence which may not be detectable given the admixture pulse scenarios fitted by Bunnefeld et al. (2018).

| No evidence for structure within Iberia
We find little variation in PSMC trajectories between different Iberian pairs in almost all species ( Figure S2) (Figure 1 and Figure S2) and so likely reflect coverage variation between individuals. E. brunniventris is the only species that showed clear variation in PSMC trajectories between Iberian pairs. While our additional analyses for E. brunniventris suggest other potential reasons for this finding (see Discussion), we cannot rule out that this variation is in part driven by population structure.

| Sensitivity analyses
Both demographic inference methods used here assume selective neutrality but make different simplifying assumptions about recombination: while PSMC approximates recombination as a Markov process along the genome, the single-step-change analyses based on the bSFS assumes no recombination within blocks. To check the extent to which recombination and selection may bias parameter estimates, we fitted histories of a step change to data simulated with recombination. Our exploration of simulated data shows that both N e and are underestimated with increasing recombination rate whereas the time of step change (T) is biased downwards when N e declines towards the present ( < 1) and upwards when N e increases ( > 1) (Table S2). These biases are an expected consequence of the shuffling of genealogical histories within blocks in the presence of recombination, which reduces the variance in bSFS configurations.
Given a recombination rate of r = 3 × 10 − 9 (Table S3), our simulation check suggests that while T, the time of the step change may be underestimated by up to a factor of two, the ability to accurately estimate , the magnitude of the step change, is little affected.

F I G U R E 3
PSMC N e trajectories for Iberian (red) and cross-population (Iberian vs Balkans) pairs (blue) of seven species of parasitoids. Population splits are visible as a divergence between the trajectories of within-population and cross-population pairs. Vertical bars indicate the split times estimated by Bunnefeld2018 using a MCL method based on the blockwise data. The split time estimated by Bunnefeld2018 for M. dorsalis is too recent (8 kya) to be visible in this figure. The time bar along the top of the figure shows (from left to right) the Holocene (H, 11.7 kya to present), the last glacial, the Eemian interglacial (E, 115-130 kya) and the penultimate interglacial

H E H E
To investigate the potential effect of selective constraint on the single-step-change analyses, we inferred single-step-change histories separately for genic and intergenic regions. Parameter estimates based on intergenic data for O. nitidulus were broadly similar to those obtained from the full data set (Table S2): while estimates of and T were slightly lower and higher, respectively, for the full data set, as would be expected as a result of selective constraint, the estimate of was little affected.

| DISCUSS ION
We analysed genome-wide sequence variation using two contrasting inference approaches to test whether and how demographic histories vary within a guild of insect parasitoids in a single Pleistocene or not. Contrary to any simple scenario of guild-wide temporal concordance in responses to Pleistocene climatic events, we find that population size changes of species in this guild differ markedly both in direction and in timescale. In fact, our single-step-change analyses reveal significant support for maximal temporal discordance; that is, each of the four species in this parasitoid guild that show evidence for a sudden change in N e has a unique time (Table 1).
Thus, our main result of temporally discordant N e change within Iberia mirrors the finding of temporal discordance by Bunnefeld et al. (2018) in the timing of divergence and admixture for this guild on a continental scale. One could argue that the signal of temporal discordance may simply be an artefact of applying an oversimple step-change model, which is unlikely to capture the subtleties of real population size change. However, the fact that we also find evidence for temporal discordance when visualizing N e change using PSMC, which makes minimal simplifying assumptions about the shape of past demographic change, suggests a genuine lack of signal for temporal concordance in this parasitoid guild across a range of spatio-temporal scales. This general finding mirrors results of other comparative studies on sets of co-distributed taxa in Europe (Ebdon et al., 2021) and the Americas (Burbrink et al., 2016;Dasmahapatra et al., 2010).
It is interesting that the only potential signal of temporal concordance we find is for four species that show the smallest change in past population sizes and which overlap in the time of peak N e ( Figure 2). While this apparent congruence may be due to a shared background demography, the signatures of which are masked in species with more drastic changes in N e , we cannot rule out alternative explanations. In particular, the fact that PSMC assumes selective neutrality is problematic given that insect genomes are more compact than mammalian genomes (Li & Graur, 1991), and so more susceptible to the effects of linked selection. Schrider et al. (2016) have shown that selective sweeps can generate troughs in PSMC trajectories while background selection has been shown to lead to erroneous inference of population growth (Ewing & Jensen, 2016).
We therefore cannot rule out the possibility that congruent peaks in N e are an artefact of selective effects which, assuming similar genome composition, mutation and recombination rates, may lead to similarly distorted PSMC inferences. However, the fact that we have reconstructed strikingly different PSMC trajectories, most of which differ markedly from the selection-induced PSMC trajectories of Schrider et al. (2016) (in that they show large declines rather than increases in N e towards the present), suggests that it is unlikely that linked selection is the main driver of the inferred N e changes. For the single-step-change analyses, we find little difference in parameter estimates when analysing all data or just intergenic regions (Table   S2), which again argues against a major effect of selection on our inferences. Furthermore, one would expect genomes with a shorter map length (physical length × recombination rate, see Table S3) to be disproportionately affected by linked selection (Mackintosh et al., 2019;Smith & Haigh, 1974). However, if anything, we observe the In general, one may view the fact that PSMC is essentially assumption-free as an advantage over model-based inference, because it enables a straightforward visualization of past N e change.
Similarly, comparing PSMC trajectories between pairs of individuals and populations allows qualitative assessment of likely periods of divergence and admixture ( Figure 3). However, the flip-side of this flexibility is that PSMC provides no obvious route for quantitatively testing (necessarily) simple demographic hypotheses across species.
Furthermore, cross-species comparisons of PSMC trajectories are problematic even when they are focused on a single clearly defined summary such as the interval of peak N e (as we have done here), because the time discretization PSMC imposes depends on the rate of coalescence, and so differs between species.

| Changes in population size coincide with late Pleistocene climatic transitions
While we emphasize that absolute time estimates need to be inter-  (Bunnefeld et al., 2018;Lohse et al., ,2010Lohse et al., , , 2012. For example, the cluster of peak N e times coincides with a large community-wide peak in the frequency of admixture between refugia inferred by Bunnefeld et al. (2018). Similarly, the large increase in N e inferred for T. auratus at the end of the Eemian interglacial period coincides with the divergence of refugial populations inferred for this species previously (Bunnefeld et al., 2018;Stone et al., 2012). It seems plausible that both events are associated with the geographic expansion of suitable oak habitat across refugial barriers during these times Petit et al., 2002), which may have triggered range expansions in both host gall wasps and their associated parasitoids.

| Population structure within Iberia and gene flow from Eastern refugia
Population structure within southern European refugia has previously been demonstrated for some species, and it has been suggested that given its complex topography Iberia should be considered a mosaic of multiple micro-refugia rather than a single entity (Feliner, 2011;Hearn et al., 2014). However, our PSMC results suggest a complete lack of population structure in six out of seven species in the parasitoid guild ( Figure S2) and imply high gene flow within Iberia.
This is perhaps unsurprising given that gall wasp-associated parasitoid wasps (and other chalcids) are able to disperse long distances even across patchy habitats and host distributions (Compton et al., 2000;Hayward & Stone, 2006

| Eurytoma brunniventris is an outlier
Eurytoma brunniventris is an outlier in our results in several ways: it is the species with the highest genetic diversity, shows signals of population structure within Iberia and is the only species for which our two inference approaches disagree. While the single-step-change analyses give no support for a change in N e , the PSMC trajectory shows a steady increase. To explore the sensitivity of the single-step-change analysis to detect gradual changes in population size, we simulated 100 replicate data sets (assuming the same size and block length as the real data) under the gradual change in N e inferred via PSMC for E. brunniventris. We find that in each case, a single-step-change history fits the data significantly better (using a parametric bootstrap analogous to that performed on the real data, see Methods) than a null model of constant N e , indicating that the single-step-change analysis is indeed sensitive to gradual increases in population size.
A possible explanation for the discrepancy between the PSMC and single-step-change analyses for E. brunniventris is that its genome assembly, the least contiguous among our set of taxa (Table S3)

| Outlook
A general question for geographically widespread communities is the extent to which component species continue to interact during the assembly process (Agosta & Klemens, 2008;Janzen, 1985;Ricklefs, 2015). Where communities are characterized by strong dependencies between species, such as specific trophic or symbiotic interactions, we might expect co-dispersal and coupled population dynamics to result in similar demographic histories (e.g. Gaume et al., 2000); this scenario is compatible with ongoing selective effects of species on each other and potential coevolution (Hall et al., 2020;Wade, 2007).
In contrast, we might expect much lower demographic concordance for members of communities characterized by weaker and less specific species interactions and only diffuse coevolution (Hall et al., 2020). This is the pattern emerging for the parasitoid assemblages attacking oak gall wasps. Though comprising a consistent set of interacting taxa spanning thousands of kilometres of longitude, the component species show highly diverse histories of range expansion (Bunnefeld et al., 2018) and changes in population size. This finding is consistent with the ability of most of the parasitoid species to attack multiple gall wasp hosts (Askew et al., 2013), weakening both direct interactions between the trophic levels and competitive interactions between parasitoid species. While there is strong evidence that host gall wasp traits structure-associated assemblages of parasitoid enemies (Bailey et al., 2009) and of parasitoid specializations for exploitation of particular hosts (Weinersmith et al., 2017), there is little evidence for specific coevolution. The population histories of very few other multitrophic communities have been explored, and the extent to which the demographic diversity within the gall wasp system is typical of parasitoid-host and other biological systems remains unknown. Work on the arthropod communities associated with the pitcher plant Sarracenia alata (Satler & Carstens, 2017) also shows variation in component species histories; further, this study suggests that more ecologically dependent associates show closer demographic and phylogeographic concordance with the host plant than less ecologically dependent species. An obvious question is whether any of the variation in patterns of N e change in our sampled parasitoid species can be attributed to variation in ecological traits (Papadopoulou & Knowles, 2016). Intriguingly, the four parasitoid species with support for a step change (in either direction) have a narrower host range (Askew et al., 2013) and a lower ancestral N e than those for which a null model of constant N e could not be rejected (Table S1). Although the number of species studied here is clearly insufficient for any statistical test of an association between host range and demographic history, these trends are compatible with the idea that ecologically specialist species (whose biology is critically dependent on interactions with a small number of other taxa) experience greater and/or more frequent changes in N e than generalists (whose demography is less tightly coupled to abundance of any specific interaction) (Östergård & Ehrlén, 2005;Rand & Tscharntke, 2007

ACK N OWLED G EM ENTS
We thank Lynsey Bunnefeld for help in the molecular laboratory and useful discussions and José Luis Nieves-Aldrey and Juli Pujade-Villar for contributing samples and Sam Ebdon for comments on the manuscript. This work was supported by a Natural Environment Research Council (NERC) grant to GNS and KL (NE/J010499/1). WW was supported by an E3 DTP studentship from NERC. KL is supported by a NERC fellowship (NE/L011522/1) and an ERC starting grant (ModelGenomLand).

AUTH O R CO NTR I B UTI O N S
KL and GS designed the project; WW analysed the sequence data with contributions from KL; all authors wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
• Raw reads have been deposited in the European Nucleotide