Comparative genetics of invasive populations of walnut aphid, Chromaphis juglandicola, and its introduced parasitoid, Trioxys pallidus, in California

Abstract Coevolution may be an important component of the sustainability of importation biological control, but how frequently introduced natural enemies coevolve with their target pests is unclear. Here we explore whether comparative population genetics of the invasive walnut aphid, Chromaphis juglandicola, and its introduced parasitoid, Trioxys pallidus, provide insights into the localized breakdown of biological control services in walnut orchards in California. We found that sampled populations of C. juglandicola exhibited higher estimates of genetic differentiation (FST) than co‐occurring populations of T. pallidus. In contrast, estimates of both the inbreeding coefficient (GIS) and contemporary gene flow were higher for T. pallidus than for C. juglandicola. We also found evidence of reciprocal outlier loci in some locations, but none showed significant signatures of selection. Synthesis and applications. Understanding the importance of coevolutionary interactions for the sustainability of biological control remains an important and understudied component of biological control research. Given the observed differences in gene flow and genetic differentiation among populations of T. pallidus and C. juglandicola, we suspect that temporary local disruption of biological control services may occur more frequently than expected while remaining stable at broader regional scales. Further research that combines genomewide single nucleotide polymorphism genotyping with measurements of phenotypic traits is needed to provide more conclusive evidence of whether the occurrence of outlier loci that display significant signatures of selection can be interpreted as evidence of the presence of a geographic mosaic of coevolution in this system.

One of the key difficulties in finding evidence for coevolution is the identification of the ecologically relevant traits or genes that are under reciprocal selection in antagonistic interactions. Recently, population genomics has been utilized to facilitate our understanding of evolutionary (Black et al. 2001, Luikart et al. 2003, Stinchcombe & Hoekstra 2008Deagle et al., 2011;Barker, Andonian, Swope, Luster, & Dlugosch, 2017) and coevolutionary processes (Parchman, Buerkle, Soria-Carrasco, & Benkman, 2016;Vermeer, Dicke, & de Jong, 2011;Yoder, 2016). Through the examination of large numbers of neutral markers, population genomics can be used to separate locus-specific effects that may be linked to genes under selection, from genomewide effects driven by genetic drift, migration, and inbreeding. This has the advantage that a population genomics approach can be applied to a wide variety of nonmodel organisms under field conditions. To investigate whether coevolutionary processes operate at different spatial or temporal scales in an antagonistic interaction, it is necessary to find evidence of the occurrence of geographic selection mosaics, trait remixing, and hot and cold spots of coevolution (Gomulkiewicz et al., 2007;Thompson, 2005). In this context, Vermeer et al. (2011) suggest that while population genomics cannot be used to test for the existence of a geographic mosaic of coevolution per se, it can be a valuable approach for the detection of unusual levels of variation or "outliers" at specific loci that are potential indicators of hot spots of reciprocal selection, and for estimation of the extent of gene flow and inbreeding that are factors contributing to trait remixing.
One field setting in which insect host and parasitoid coevolution is thought to play an important role (see Holt & Hochberg, 1997) is during the importation and establishment of non-native parasitoids to suppress the abundance of invasive insect pests (Heimpel & Mills, 2017;Hoddle, 2004;Van Driesche et al., 2010). Biological control programs are known to be well suited for the study of evolution (Roderick, Hufbauer, & Navajas, 2012;Roderick & Navajas, 2003), but few studies exist that have used population genetic techniques to conduct comparative analyses of evolutionary change postintroduction. This may, in part, be explained by a predominant focus in biological control on pre-introduction surveys for natural enemies without sufficient emphasis on longer-term postintroduction monitoring (McCoy & Frank, 2010;Mills, 2000Mills, , 2017).
Our study system consists of walnut, an exotic tree crop in California; walnut aphid, Chromaphis juglandicola (Kaltenbach), an invasive species that is active from March until early December (Sluss, 1967); and Trioxys pallidus (Haliday), an introduced exotic parasitoid wasp. As is typical for nonhost alternating aphids, walnut aphids reproduce through cyclical parthenogenesis (Simon, Rispe, & Sunnucks, 2002), in which females produce multiple generations of female offspring asexually through the summer. In the fall, decline in photoperiod triggers the development of a single sexual generation of both male and female aphids, and oviparous females deposit eggs which overwinter until the following spring (Davidson, 1914). Walnut aphid appears not to have secondary defensive symbionts (Russell, Latorre, Sabater-Muñoz, Moya, & Moran, 2003), and our own surveys of populations in California support this earlier observation (J.C. Andersen, unpublished data). Biparental hymenopteran parasitoids, such as T. pallidus, reproduce through haplodiploidy in which haploid males have a single maternally inherited copy of each chromosome and diploid females have both a paternal and maternal copy of each chromosome (Heimpel & de Boer, 2008). This wasp was originally introduced from southern France, which resulted in establishment in the southern and coastal regions, but failed to establish the parasitoid in the primary walnut growing region of the Central Valley (Schlinger, Hagen, & van den Bosch, 1960).
Subsequently, a second introduction from Iran led to widespread establishment and reduction in walnut aphid densities throughout California (van den Bosch et al., 1979). While outbreaks of walnut aphids have occurred since the establishment of the Iranian strain of T. pallidus, these were associated with the use of azinphosmethyl, an insecticide used for the control of other walnut pests, and a resistant population of T. pallidus was reared and released (Brown, Cave, & Hoy, 1992;Hoy & Cave, 1989;Hoy et al., 1990). More recently, localized increases in the abundance of aphids have led to a resumption of in-season insecticidal treatments in walnut orchards in California (Hougardy & Mills, 2008), and the reason for this remains unknown.
In a previous study, we examined whether hybridization among descendants of two different introduced populations of T. pallidus may have played a role in the observed breakdown of biological control services (Andersen & Mills, 2016). While hybridization was found to be rare in California, we did find evidence of genetic differentiation among populations of the introduced parasitoid. Given the genetic structuring of T. pallidus populations in walnut orchards in California, we were interested to know whether C. juglandicola populations displayed similar patterns of differentiation and how these patterns varied geographically within California. Therefore, the objectives of this study were (1) to compare levels of genetic differentiation among populations of T. pallidus and C. juglandicola in California, (2) to estimate rates of gene flow and inbreeding among these populations for both species, and (3) to use population genetics to detect outlier loci and the potential for reciprocal selection as preliminary evidence for the existence of a geographic mosaic of coevolution in the walnut aphid biological control program.

| Sampling locations
Californian walnut orchards were visited between 2010 and 2014.
Orchards were selected to represent a broad range of geographic locations, but there was no prior information on the history of walnut aphid densities or levels of parasitism in each orchard. At each location, we collected individuals identified as T. pallidus either by aspirating adults or by collecting mummified walnut aphids and placing small cut-out sections of leaf material with each mummy into glass vials (9.5 mm × 3 mm).
These vials were closed with a foam stopper and stored at room temperature until adults emerged. Whether aspirated, or reared, adults of T. pallidus were then stored in 95% ethanol at −20°C for molecular analysis.
Individuals of C. juglandicola were collected in the field and immediately placed in 95% ethanol and then stored at −20°C for molecular analysis.
Effort was taken to collect only a few individuals per tree and to prioritize sampling from as many different trees as possible in each orchard to reduce the sampling of clonally related individuals (Lozier, Roderick, & Mills, 2007). Full details for each of the sampling locations are presented in Table 1.

| DNA extraction and microsatellite genotyping
DNA was extracted from adult females of T. pallidus and C. juglandicola by grinding individuals with a mortar and pestle and then using the modified DNA extraction protocols presented in Andersen and Mills (2014). Standard PCR protocols were then used to amplify 15 polymorphic microsatellite markers for T. pallidus and 12 polymorphic microsatellite markers for C. juglandicola following protocols presented in Andersen and Mills (2014). Briefly, microsatellite loci were amplified from 7 to 20 aphids and/or parasitoid females from each location using fluorescently labeled primers, and PCR products for up to four loci were pooled before genotyping so that no two loci with the same fluorescent label were combined. Products were then genotyped on an Applied Biosystems 3730XL DNA Analyzer at the University of California Berkeley DNA Sequencing Facility using the LIZ 600 size standard, and fragment lengths were then scored using the Microsatellite Plug-in for Geneious Pro v. 5.6.2 (Drummond et al., 2012). the presence of locus-by-locus linkage disequilibrium (LD) were then estimated with the software package GenePop (Raymond & Rousset, 1995;Rousset, 2008). Estimates of population differentiation based on F ST were generated using FreeNA (Chapuis & Estoup, 2007) to account for the potential presence of null-alleles, and whether populations were significantly differentiated between each population pair was determined using the exact G test implemented in GenePop. Recent migration rates (i.e., the proportion of individuals in a population that were estimated to be derived from a second population)

T A B L E 1 Collection locality and genetic summary information for populations of Trioxys pallidus and Chromaphis juglandicola
were then estimated between each population using the BayesAss Edition v. 3.0 (BA3) software package (Wilson & Rannala, 2003). Four independent analyses for each species were conducted, each using a mixing parameter of 0.8 for allele frequencies, migration rates, and inbreeding coefficients and a runtime of 10 million generations with a burn-in period of 1 million generations. Results were then visualized and summarized across runs for each species using the program Tracer v. 1.6.0 (Rambaut & Drummond, 2007). Finally, evidence for population structure (genetic diversity) and outlier loci can be used to identify potential hot and cold spot locations, which in conjunction with measurements of phenotypic or behavioral traits, and can be used to confirm whether the interacting species are under reciprocal selection. By combining measures of gene flow and inbreeding with outlier detection using both known cold spot and undetermined hot and cold spot locations, this approach addresses the three underlying processes that form the basis for testing the geographic mosaic theory of coevolution; (1) coevolutionary hot and cold spots, (2) selection mosaics, and (3) trait remixing (Gomulkiewicz et al., 2007).

| Geographic mosaic of coevolution
Following this approach, we utilized two of our sampled orchards from which only one of the two interacting species was present to act as our reference coevolutionary cold spot locations.  (Luu, Bazin, & Blum, 2017), it has been widely used for the analysis of microsatellite datasets.

| Population genetic analyses
Genotyping results for the individuals analyzed are available in Appendix S1.

| Detection of outlier loci
Analyses based on identifying outliers from values of null-allele corrected F ST for each population pair using the "boxplot" function in R identified outliers in three of the pairwise comparisons for T. pallidus Using the approach of Vermeer et al. (2011), results suggested that Escalon, Newark, and Upper Lake potentially represent coevolutionary hot spots, while Arbuckle, Linden, and Yuba City potentially represent coevolutionary cold spots. However, Bayesian simulations using BayeScan indicated that none of the potential outliers for either species showed significant (p < .05) signatures of selection.
As such, coevolution has long been thought to play an important role in the sustainability of biological control services (Holt & Hochberg, 1997;Jones, Vanhanen, Pettola, & Drummond, 2014;Kraaijeveld & Godfray, 1999 [Davidson, 1914;van den Bosch et al., 1979]). Alternatively, it could also be due to differences in the rates of evolution as a result of selection and/or to differences in reproductive strategies (i.e., sexual and haplodiploid for T. pallidus versus cyclical parthenogenesis and diploid for C. juglandicola). While aphid species have shown evidence of rapid evolution in response to changes in their environment (e.g., Harmon, Moran, & Ives, 2009), a recent study of aphid-parasitoid coevolution found evidence for genetic tracking of both species (Nyabuga et al., 2012). In this latter study, the authors also found greater levels of differentiation among aphid populations compared to their parasitoids. However, in contrast to our results, they found that the aphid populations had greater inbreeding coefficients than the parasitoid populations and suspected that this arose from different metapopulation dynamics. The authors also considered that the rate of evolution of the parasitoid relative to that of its aphid host was constrained by the lag time in colonization of new patches (Nyabuga et al., 2012). Lag time may disrupt reciprocal selection (Lapchin & Guillemaud, 2005), and a difference in evolutionary rates of interacting organisms can have negative impacts on the stability of their relationships. For example, predator-prey dynamics can be negatively affected by the rapid evolution of the prey species (Yoshida, Jones, Ellner, Fussmann, & Hairston, 2003), while conversely, herbivore-plant dynamics can be negatively affected by the rapid evolution of the herbivore (Smith, de Lillo, & Amrine, 2010). However, why T. pallidus with its greater levels of gene flow among populations would display higher levels of G IS is unclear as mathematical models predict that as migration rates increase among populations, local adaptation within those populations will decrease (Blanquart, Gandon, & Nuismer, 2012).
In a pioneering paper, Holt and Hochberg (1997)  documented or is suspected to have occurred, in at least a couple of biological control programs (e.g., Goldson et al., 2014;Ives & Muldrew, 1984;Tomasetto, Tylianakis, Reale, Wratten, & Goldson, 2017). In one of these programs, the control of the larch sawfly in Canada, resistance to parasitism may have arisen due to the accidental importation and spread of pre-adapted resistant host strains (Ives & Muldrew, 1984). For another, the control of the Argentine stem weevil in New Zealand, parasitism rates declined by nearly 50% over a 5-year period (Goldson et al., 2014) starting exactly 7 years after parasitoid release irrespective of the actual year of introduction (Tomasetto et al., 2017).
Similar to the walnut aphid biological control program, the Argentine stem weevil and its introduced parasitoid differ in their reproductive strategies (sexual for the Argentine stem weevil versus asexual for the parasitoid). Therefore, it is possible that coevolutionary interactions in biological control systems that rely on natural enemies with a reproductive strategy that differs from that of their target host may become decoupled due to different rates of evolution among the interacting species.

| Geographic mosaic of coevolution
If as expected, coevolutionary interactions are both spatially and temporally dynamic (Torres, 2009) and are important for biological control services (Holt & Hochberg, 1997;Jones et al., 2014;Kraaijeveld & Godfray, 1999, then localized breakdowns in biological control services might be a common and transitory occurrence as predicted by the geographic mosaic theory of coevolution (Thompson, 1994(Thompson, , 2005. While geographic mosaics have been observed for interactions between herbivores and plants (Muola et al., 2010;Siepielski & Benkman, 2005;Vermeer, Verbaarschot, & de Jong, 2012), predators and prey (Brodie & Ridenhour, 2002), and hosts and parasites (Dixon, Craig, & Itami, 2009;Lorenzi & Thompson, 2011;Thompson, 2009b This is in part due to the fact that there have been very few longterm, postrelease studies of biological control agents and their targets (McCoy & Frank, 2010;Mills, 2000Mills, , 2017, and that it can be difficult to identify a priori which adaptive traits to measure when studying coevolutionary interactions. In this context, comparative population genomics may provide a useful approach to obtain preliminary evidence for the presence and/or potential importance of coevolution in biological control systems. Based on our findings, we suspect that coevolution is important for the sustainability of biological control programs and that long-term studies would likely reveal a continuum from sustained effective control when coevolutionary interactions are strong, to failures when they are weak. Under this scenario, biological control programs may experience temporary failures in effective control at a localized scale, and yet experience sustainable control at a regional or landscape scale due to connectivity and movement between local populations as predicted by the geographic mosaic theory of coevolution.

ACKNOWLEDGMENTS
Funding for this study was made possible by an EPA-STAR Fellowship and a Robert van den Bosch Memorial Fellowship awarded to JCA.
The authors would like to thank K. Anderson, and M. Labbé for their field assistance in collecting T. pallidus and C. juglandicola. In addition, we would like to thank the University of California Cooperative Extension, and the East Bay Regional Parks Department, as well as the individual walnut growers who granted us access to their properties.
We would also like to thank J. Lozier and B. Ort for their technical assistance in conducting the molecular analyses; and R. Nielsen, G.
Roderick, C. Smith, and two anonymous reviewers for their comments on earlier drafts of this manuscript.

CONFLICT OF INTEREST
None declared.

AUTHORS' CONTRIBUTIONS
JCA conduct the field sampling, generated the molecular data, and performed analyses. NJM oversaw the collection of data and provided laboratory support. Both authors contributed equally to the study design and the drafting of the manuscript.