Habitat or temporal isolation: Unraveling herbivore–parasitoid speciation patterns using double digest RADseq

Abstract Ecological speciation is often observed in phytophagous insects and their parasitoids due to divergent selection caused by host‐associated or temporal differences. Most previous studies have utilized limited genetic markers or distantly related species to look for reproductive barriers of speciation. In our study, we focus on closely related species of Lygus bugs and two sister species of Peristenus parasitoid wasps. Using mitochondrial DNA COI and genomewide SNPs generated using ddRADseq, we tested for potential effects of host‐associated differentiation (HAD) or temporal isolation in this system. While three species of Lygus are clearly delineated with both COI and SNPs, no evidence of HAD or temporal differentiation was detected. Two Peristenus sister species were supported by both sets of markers and separated temporally, with P. mellipes emerging early in June and attacking the first generation of Lygus, while P. howardi emerging later in August and attacking the second generation of their hosts. This is one of the few studies to examine closely related hosts and parasitoids to examine drivers of diversification. Given the results of this study, the Lygus‐Peristenus system demonstrates temporal isolation as a potential barrier to reproductive isolation for parasitoids, which could indicate higher parasitoid diversity in regions of multivoltine hosts. This study also demonstrates that incorporating systematics improves studies of parasitoid speciation, particularly by obtaining accurate host records through rearing, carefully delimiting cryptic species and examining population‐level differences with genomic‐scale data among closely related taxa.


| INTRODUC TI ON
A growing number of evolutionary studies have focused on ecological speciation, in which new species arise as a result of ecologically driven divergent selection Hood et al., 2015;Nosil, Crespi, & Sandoval, 2002;Rundle & Nosil, 2005;Schluter, 2009). During ecological speciation, reproductive barriers arise as a by-product of adaptation to divergent environments. Ecological speciation has been observed in herbivorous insects in the form of host-associated differentiation (HAD), where specialists diverge through phenological or host shifts as a result of competition and/ or predation (Nosil et al., 2002;Rundle & Nosil, 2005;Schluter, 2009). Adaptation to divergent host plants leads to an accumulation of multiple reproductive barriers, ultimately resulting in the separation and formation of new species (Dres & Mallet, 2002;Forbes et al., 2017). The presence of HAD is often associated with a temporal component, where temporal divergence in the breeding time over timescales ranging from days, seasons, or even years (Taylor & Friesen, 2017). Temporal isolation can contribute to divergence alone or concurrently with traits such as host preference to reinforce divergence along the speciation continuum Feder et al., 1994;Taylor & Friesen, 2017). Although allopatric populations are often defined by spatial differentiation, populations with overlapping distributions and phenological differences can also be argued as allopatric, but on a temporal scale (Taylor & Friesen, 2017). Most documented cases of temporal speciation among phytophagous insects involve seasonal separation of breeding time after host shifts resulting in selection of synchrony with host phenology, contributing to reproductive isolation as selection Feder et al., 1994;Nosil et al., 2002;Stireman, Nason, & Heard, 2005). These phenological shifts are often associated with genes controlling diapause duration, timing of diapause termination, and circadian rhythms, which could contribute to divergent selection that ultimately drives ecological speciation (Ragland, Egan, Feder, Berlocher, & Hahn, 2011;Ragland, Sim, Goudarzi, Feder, & Hahn, 2012;Ragland et al., 2017;Taylor & Friesen, 2017).
Numerous studies have shown that generalist insect herbivore "species" are often multiple, genetically divergent cryptic lineages, each specializing on a subset of the full host plant range (Dres & Mallet, 2002;Peccoud, Ollivier, Plantegenest, & Simon, 2009;Powell, Forbes, Hood, & Feder, 2014). This is an important distinction as true generalists feed on a variety of host plants indiscriminately, while cryptic specialists exhibit host preferences that were overlooked due to morphological similarities. Therefore, the accurate identification of true generalists from cryptic specialists in various stages of speciation is vital to studies on the effects of host or temporal differences on biogenesis.
In addition, HAD can have rippling effects at higher trophic levels, resulting in divergence of parasitoids in the form of cascading/sequential HAD (Abrahamson & Weis, 1997;Forbes, Powell, Stelinski, Smith, & Feder, 2009;Hood et al., 2015;Nicholls, Schönrogge, Preuss, & Stone, 2018;Stireman, Nason, Heard, & Seehawer, 2006). As many parasitoids are also cryptic specialists that are tightly linked to the phenology of their hosts, cascading HAD on species lineages of herbivores could result in the sequential radiation of these hyperdiverse lineages of parasitoids (Forbes et al., 2009;Hood et al., 2015;Stireman et al., 2006). However, many previous studies of HAD and sequential HAD were limited to few molecular markers (Antwi et al., 2015;Hood et al., 2015;Leppanen et al., 2014;Nicholls et al., 2018;Stireman et al., 2006), which provides limited molecular characters to examine fine-scaled species-level differentiation.
In addition, most studies focus on specialist herbivores with few studies on parasitoids. Studies that have involved examinations of parasitoids have mainly included assemblages of distantly related parasitoids that make inferences about drivers of diversification in upper trophic levels difficult Nicholls et al., 2018;Stireman et al., 2006). Therefore, studies focusing on closely related parasitoids species are needed to examine patterns of speciation due to ecologically divergent selection.
Accurate delimitation of divergent lineages is paramount to speciation studies, as they are often morphologically cryptic. Studies utilizing variations in restriction-site associated DNA sequencing (RADseq) to delimit species and determine drivers of divergence have become more abundant (Bagley, Sousa, Niemiller, & Linnen, 2017;Bernal, Gaither, Simison, & Rocha, 2017;de Oca et al., 2017;Eaton & Ree, 2013). RADseq approaches are less susceptible to incomplete lineage sorting and introgression than traditional multigene methods (Andrews, Good, Miller, Luikart, & Hohenlohe, 2016). This method is ideal for detecting population/species-level differences and has been shown to be promising for studies on ecological speciation of herbivorous insects (Bagley et al., 2017;Egan et al., 2015).
Studying the reproductive barriers of parasitoid in relation to their hosts is central to understanding origins of parasitoid diversity and may also provide important insights into conservation biology as parasitoids have been shown to provide ecosystems with trophic redundancy that reduces extinction risks (Sanders, Thébault, Kehoe, & van Veen, 2018). In addition, understanding the intimate relationships between pestiferous herbivores and their parasitoids would greatly improve the success rate of biological control programs (Peixoto et al., 2018;Zhang, Ridenbaugh, & Sharanowski, 2017). To that end, we investigate potential reproductive barriers in the Lygus-Peristenus system, which includes a genus of economically important herbivores and the parasitoid species that attack them.
The herbivores in this system are plant bugs in the genus Lygus Hahn (Hemiptera: Miridae), which include many species of generalist agricultural pests (such as Lygus lineolaris Palisot de Beauvois) that feed on a variety of economically important crops. Although HAD has been recorded from other Miridae (Hereward, Walter, Debarro, Lowe, & Riginos, 2013), no evidence of HAD has been shown in Lygus species despite the detection of population-level differences based on geography (Burange, Roehrdanz, & Boetel, 2012;Zhou, Kandemir, Walsh, Zalom, & Lavine, 2012). Lygus have one to three generations per year depending on the temperature, where southern populations in warmer climates are multivoltine and northern populations in cooler climates tend to be univoltine (Cárcamo et al., 2002;Haye et al., 2013). The Canadian prairies ecosystem is a major agricultural growing region where Lygus is an economically relevant pest on several field crops, such as canola, alfalfa, and mustard.
Closely related species are often found in sympatry, so HAD may be a driver of population divergence in this system, as populations could be cryptic, specializing on certain plants.
Species of Peristenus (Hymenoptera: Braconidae) are widely distributed koinobiont endoparasitoids of nymphal plant bugs, including Lygus species (Zhang, Stigenberg, Meyer, & Sharanowski, 2018). A recent revision of the Nearctic Peristenus pallipes complex synonymized nine species recognized by Goulet and Mason (2006) to just three based on morphometrics, mitochondrial DNA (COI and CytB), and ecological differences (Zhang et al., 2017). This revision also demonstrated a range overlap for Peristenus dayi Goulet with sister species Peristenus mellipes (Cresson) and Peristenus howardi Shaw in southern Alberta (Zhang et al., 2017). As these Peristenus species persist in sympatry, there are likely reproductive barriers preventing hybridization and interbreeding between species. These may be ecological isolating mechanisms, such as differences in micro-habitat, emergence timing, and reproduction. Peristenus host preference may also explain the maintenance of three sympatric species, but due to morphological similarity among Lygus nymphs, host records are often listed simply as Lygus species (Goulet & Mason, 2006). The drivers and maintenance of species boundaries in these closely related parasitoids are unknown, but a likely explanation is divergence through sequential HAD as their hosts specialize and diverge.
In this study, we used a combination of COI (mtDNA) and double digest RADseq (ddRADseq) (Peterson, Weber, Kay, Fisher, & Hoekstra, 2012) to test for barriers of reproductive isolation in closely related parasitoids. We (a) confirm monophyly and delimit species of Lygus and their Peristenus parasitoids; (b) test for potential host plant associations or temporal differentiation on sympatric species of Lygus; and (c) determine whether sequential HAD or temporal differentiation are driving forces of speciation on sympatric species of Peristenus. As herbivore-parasitoid evolutionary histories can provide valuable insights into the genesis of biodiversity, this is one of the first studies to address the evolutionary patterns within a tritrophic system that utilizes host plant, herbivore, and parasitoid using next-generation sequencing data and closely related parasitoids.

| Sample collection and DNA Extraction
To obtain Peristenus with accurate host records delineated to species, we sampled early instar nymphal Lygus bugs weekly from May to August of 2015 from two sites in Lethbridge, Alberta, as this is the only region in which the range of both P. mellipes and P. howardi overlaps (Sharanowski, Zhang, & Wanigasekara, 2014;Zhang et al., 2017). One additional site where only P. mellipes is found was sam- If the Lygus nymphs were parasitized, the emerged larval parasitoid and dead host were preserved in 95% EtOH until DNA extraction. Genomic DNA was extracted following the DNeasy Tissue Kit Protocol (Qiagen, Valencia, CA, USA), using a destructive sampling method as the larval parasitoid and host nymphs were unidentifiable morphologically. We quantified the concentration of DNA extracts using Quant-iT High-Sensitivity DNA Assay Kit (Invitrogen, Eugene, OR, USA). Peristenus dayi was excluded from this study despite being closely related to the other parasitoids, as it parasitizes Adelphocoris lineolatus (Goeze), a distant relative of Lygus within Miridae, and we were interested in patterns between closely related herbivores and parasitoids.
We used a modified ddRADseq protocol from Peterson et al.  We used ipyrad v0.7.23 (Eaton, 2014) to process raw sequences, using the following stringent settings to ensure the data quality for downstream analyses after parsing out Lygus from Peristenus:

| Phylogenetic analyses
The best-fitting model of molecular evolution for COI was tested using A maximum-likelihood supermatrix approach using the concatenated ddRADseq SNPs dataset was also conducted with RAxML 8.2.0 (Stamatakis, 2006), using the GTR + Γ model of nucleotide substitution and 1,000 bootstrap pseudoreplicates. The resulting trees were visualized and modified in the same manner as the COI trees.

| Population genomic analyses
To determine whether there was population structure within clades identified in the phylogenetic analysis, we performed a Bayesian clustering analysis for both Lygus and Peristenus unlinked SNP datasets (1 SNP per locus) from the ipyrad output stated earlier without prior assignments in Structure v 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). Ten runs were completed for each population (K) up to the maximum number of populations within each clade using 100,000 burn-ins and 500,000 replicates for each run. The R package pophelper (Francis, 2017) was used to visualize the diagrams.
The Evanno ΔK method (Evanno, Regnaut, & Goudet, 2005) was used in Structure Harvester v 0.6.94 (Earl, 2012) to determine the most likely value for K. We also created a custom dataset of the SNPs containing only Alberta populations of P. mellipes and P. howardi in ipyrad using the same settings discussed above. We tested for potential genetic differences under selection between the Alberta populations where the two Peristenus species are found in sympatry.
Impacts of locality, host association, and time of emergence on genetic variation of the three Lygus species were tested using AMOVA (analysis of molecular variance) using clustering between localities (for L. borealis), host plants (for L. keltoni and L. elisus), and collecting dates for all three species of Lygus. Similarly, AMOVA was used to test for differences between hosts for both Peristenus species and difference between collection localities for P. mellipes. All AMOVAs were conducted with R packages adegenet (Jombart & Ahmed, 2011) and poppr (Kamvar, Tabima, & Grünwald, 2014) using the full SNP dataset as described above.  Figure S1) and both species of Peristenus (Supporting information Figure S2).

| Phylogenetic analyses
A total of 33 Lygus (543 bp) and 37 Peristenus (629 bp) COI sequences were used for the phylogenetic analyses (Table 1). Three monophyletic clades of Lygus were identified based on the monophyletic clustering with identified specimens available in BOLD: Lygus borealis (Kelton), Lygus keltoni Schwartz, and Lygus elisus Van Duzee (Supporting information Figure S1). All three species of Lygus were collected in Alberta, while only L. borealis was collected in Manitoba.
Both L. keltoni and L. elisus were collected from all three host plants, while L. borealis was collected exclusively on alfalfa (Table 1)

. Both
Peristenus mellipes and P. howardi were recovered as monophyletic clades (Supporting information Figure S2). Peristenus mellipes was reared from all three Lygus species and found in both Manitoba and Alberta, while P. howardi was reared from L. borealis and L. keltoni and was found exclusively in Alberta (Table 1; Supporting information Figure S2).

| Population genomic analyses
Using the ΔK approach, Bayesian clustering analyses in STRUCTURE indicated K = 3 (Figure 3a) in Lygus, which corresponds to the number of species identified by phylogenetic methods (Figure 2). The STRUCTURE results show K = 3 among the two Peristenus species, as population structure was not found within P. howardi, but splits P. mellipes into an Alberta-specific population and a Manitoba population ( Figure 3B).
No significant genetic differentiation was detected among any of the AMOVA partitions (locality, host plant, collecting date) for the three Lygus species (Table 2). No differences between host bugs were detected for both species of Peristenus (Table 3a), but significant genetic differences (p = 0.01) were detected among collection localities within P. mellipes, explaining 11.77% of the genetic variation (Table 3b).

| Identification of Lygus and Peristenus species using molecular data
The accurate identification of Lygus species has been problematic in the past, because of the inconsistency between morphological differences of nymphs and COI data (Gwiazdowski, Foottit, Maw, & Hebert, 2015). The Lygus species included in this study, L. borealis, L. elisus, and L. keltoni, were often misidentified even by experts because of their variable adult phenotypes (Gwiazdowski et al., 2015). This taxonomic confusion has made previous host plant records in this group TA B L E 2 Analysis of molecular variance (AMOVA) using clustering between (a) localities, (b) host plants, and (c) collecting dates for all three species of Lygus used in this study

| Lack of HAD and Temporal isolation within Lygus species
Based on our phylogenetic analyses on Lygus ( Figure 1) and AMOVA (Table 2), it is unlikely that Lygus species evolved through hostassociated differentiation in the Canadian prairies. The three species of Lygus are all generalist herbivores feeding on a variety of available food sources, as no host plant-specific lineages were found within each species ( Figure 1, Table 2). While both L. elisus and L. keltoni were found on all three host plants sampled in this study, L. borealis were only found from alfalfa. The apparently narrow host range of L. borealis could be a by-product of our sampling, as they have been collected from other host plants such as canola (Brassica spp.) in other studies (Cárcamo et al., 2002;Otani & Cárcamo, 2011). These results show that Lygus species are truly generalists as we found no genetic divergence based on host. This lack of HAD is consistent with studies of other Lygus species such as L. lineolaris (Burange et al., 2012) and L. hesperus (Zhou et al., 2012) despite the detection of population-level differences, indicating that factors other than HAD likely drove their evolution.

| Temporal isolation but no HAD within the Peristenus species
Peristenus host choice was not significantly different in the hierarchical AMOVAs (  (Zhang et al., 2017). Our findings are consistent with Fernández, Laird, Herle, Goulet, and Cárcamo (2018), who found P. mellipes occurs early in the season between late May and late July and P. howardi in late June to late August. In addition, emergence times of P. mellipes were on average 13 days earlier than P. howardi in laboratory trials (Fernández et al., 2018). It is unknown how frequently parasitoids exhibit temporal speciation, as most of previous works on ecological speciation have focused on herbivorous insects (Forbes et al., 2017). However, the development of reproductive isolation as a by-product of divergent ecological selection should have similar genomewide effects as herbivorous insects, especially if considerable standing genomic variation is already present Michel et al., 2010).
However, most of the genetic variation is still within samples of each site (55.78% variation, p = 0.01), suggesting that other factors are responsible for the genetic variation observed. Additionally, no hostassociated patterns were observed as Manitoba samples only consisted of wasps reared from L. borealis feeding on alfalfa (Table 3).
The Manitoba P. mellipes has only one generation per year despite the absence of P. howardi, which could be the result of their host phenology as Manitoba has a shorter summer than Alberta, thus only allowing for the development of one full generation of Lygus (Haye et al., 2013). While P. mellipes were only collected from Canadian prairies in this study, previous work (Zhang et al., 2017) and historical records have shown that there are two generations of Lygus and P. mellipes in warmer regions such as Ontario (Goulet & Mason, 2006

| CON CLUS ION
Using mitochondrial DNA and genomewide SNPs, our comparative analysis of genetic differentiation between the two sister Peristenus species attacking multiple Lygus hosts revealed temporal divergence rather than host-associated differentiation. Temporal isolation likely played a vital role in the speciation process of Peristenus, whether it is acting alone or in concert with host preferences or other preor postzygotic barriers to gene flow. This is one of the first studies to demonstrate the potential of genomic data in resolving the tritrophic evolutionary relationships between plant, herbivore, and parasitoids. This study also demonstrates the importance of systematics to studies of parasitoid speciation, particularly careful delimitation of cryptic species, host rearing to obtain accurate records, and genomic-scale data for examining any population-level differences among closely related taxa.
Given the results of our study, the Lygus-Peristenus system can also be added to the growing body of literature on the importance of temporal separation as a driving force for ecological speciation and its effect on the evolution of the rich diversity of life. Currently, the importance of temporal differences in parasitoid speciation is poorly understood, but temporal isolation likely plays a significant role in the adaptations to host phenology.
Many phytophagous insects and their parasitoid systems are well studied because of their agricultural and economical importance; thus, large, collaborative, genomic-scale studies exploring these taxa could yield valuable insights into the prevalence and impact of temporal isolation in host-driven ecological speciation.

ACK N OWLED G EM ENTS
We would like to thank Leah Irwin and Melanie Scallion for assistance with specimen collecting and parasitoid rearing; Patrick Larabee for debugging and computational issues; and Michelle Gaither, Anna Forsman, and Alexa Trujillo for troubleshooting ddRADseq issues and data interpretation. We would also like to thank Robin Bagley for providing valuable feedback on earlier drafts of this paper. This research was supported in part from funding through the Western Grains Research Foundation, Growing Forward 2 program from Agriculture and Agri-Food Canada, an NSERC discovery grant and internal funding from the University of Central Florida, all awarded to B.J. Sharanowski.

CO N FLI C T O F I NTE R E S T
None declared.

DATA ACCE SS I B I LIT Y
COI sequences are available on GenBank Accession Numbers