Fitness and host use remain stable in a biological control agent after many years of hybridization

decrease fitness, with consequences for control efficacy. Additionally, hybridization may lead to changes in host use patterns that may put non-target host species at increased risk. We collected genomic data to determine ancestry of individuals from ten collection sites across a hybrid zone of Diorhabda spp., the biological control agent of Tamarix spp. in North America. We paired genomic data with phenotypic data on fitness proxies of body size and early fecundity measured in a common garden, and host use for individuals from three of those sites. We found two originally released pure species and a wide range of introgression in hybrids between three Diorhabda species. Body size and early fecundity were similar in pure species and hybrids, indicating hybridization is not detrimental to insect fitness or the biocontrol program and may provide variation that allows populations to become locally adapted. Host use of hybrids was very similar to that of pure species, although some hybrid individuals had increased preference for Frankenia salina , a native non-target species. We find that hybridization has likely not been detrimental to the efficacy and safety of the Diorhabda biocontrol program, but possible impacts on F. salina should be monitored, considering ongoing hybridization and evolution in the field.


H I G H L I G H T S
• We examine fitness and host use of hybrid biological control agents from the field.
• Phenotypes paired with genomics can assess the risk of non-target attack by hybrids.
• Fitness and host use are stable in hybrid and non-hybrid Diorhabda biocontrol agents.
• Variation associated with non-target host use may be higher in hybrids.

Keywords:
Introgression Hybrid breakdown Hybrid vigor Non-target effects Diorhabda Risk assessment A B S T R A C T Releasing several ecotypes or species of biological control agent is common in many biological control programs. However, the consequences of hybridization for fitness and host use of the resulting hybrids are difficult to predict, especially for hybrids between more than two species and with varying levels of introgression. Hybridization of biocontrol agents may increase or decrease fitness, with consequences for control efficacy. Additionally, hybridization may lead to changes in host use patterns that may put non-target host species at increased risk. We collected genomic data to determine ancestry of individuals from ten collection sites across a hybrid zone of Diorhabda spp., the biological control agent of Tamarix spp. in North America. We paired genomic data with phenotypic data on fitness proxies of body size and early fecundity measured in a common garden, and host use for individuals from three of those sites. We found two originally released pure species and a wide range of introgression in hybrids between three Diorhabda species. Body size and early fecundity were similar in pure species and hybrids, indicating hybridization is not detrimental to insect fitness or the biocontrol program and may provide variation that allows populations to become locally adapted. Host use of hybrids was very similar to that of pure species, although some hybrid individuals had increased preference for Frankenia salina, a native non-target species. We find that hybridization has likely not been detrimental to the efficacy and safety of the Diorhabda biocontrol program, but possible impacts on F. salina should be monitored, considering ongoing hybridization and evolution in the field.

Introduction
Hybridization can drive both ecological and evolutionary change. For example, it can increase invasiveness (Hovick and Whitney, 2014;Mesgaran et al., 2016), alter which environments a population can use (Rieseberg et al., 2003), and influence interactions between species through changes in traits such as plant secondary chemicals (Orians, 2000;Whitney et al., 2006). Given these kinds of effects, hybridization may have important consequences for classical biological control (hereafter, biocontrol), in which natural enemies are intentionally introduced for the control of invasive pest species (Heimpel and Mills, 2017;Szűcs et al., 2021;Wright and Bennett, 2018). In classical biocontrol programs, multiple ecotypes (Lowry, 2012) or even distinct species may be released (Van Driesche and Bellows, 1996) and hybridization among them may affect fitness, which influences the efficacy of control, and host use, which affects the safety of non-target species.
Hybridization may directly improve fitness through masking of deleterious recessive alleles in heterozygotes (Drake, 2006). Additionally, hybridization can facilitate adaptation to the environment by increasing phenotypic and genetic diversity (Ellstrand and Schierenbeck, 2000;Hovick and Whitney, 2014;Schierenbeck and Ellstrand, 2009). Hybridization can also reduce fitness if it disrupts local adaptation or creates deleterious gene combinations (i.e., incompatibilities). A challenge in understanding such outbreeding depression (reduction in fitness after hybridization) is that it may take several generations to manifest (Edmands, 2007(Edmands, , 2002. The effects of hybridization on fitness may also change over time, for instance fitness may increase in early generations through hybrid vigor, then subsequently decrease for several generations due to outbreeding depression, and then increase again as selection acts on the increased genetic variation. The fitness outcomes of hybridization are crucial to biocontrol agent efficacy, as populations of highly fit agents can grow rapidly, enabling them to better control the target pest. The potential effects of hybridization on host use are complex (Wright & Bennett 2018). Host use of highly specialized phytophagous insects, such as those used in the biocontrol of weeds, seems to be quite canalized and not prone to evolving rapidly or forming extreme phenotypes (i.e., transgressive segregation) (Hardy et al. 2020). Indeed, there are no examples of evolutionary shifts in the fundamental host range of biological control agents of invasive weeds, where the fundamental host range is defined as all possible hosts that can be used by an agent (Hinz et al., 2019;van Klinken and Edwards, 2002;Wright and Bennett, 2018). However, shifts in preference for or performance on hosts already within the fundamental host range have occurred in a few agents (Bitume et al., 2017;Fukano et al., 2016;Hoffmann et al., 2002;Thomas et al., 2010), often in conjunction with hybridization. For example, two biotypes of the scale insect Dactylopius tomentosus have distinct host preferences for target cactus species in the genus Cylindropuntia, but their hybrids have higher fitness than parental biotypes and no longer exhibit differences in host preference (Mathenge et al. 2010). In another example, the effects of hybridization between two biotypes of another scale insect, Dactylopius opuntiae, which prefer distinct Opuntia cactus species, depended on generation: F1 hybrids could develop on both Opuntia species, but half of F2 hybrids reverted to being host specific (Hoffmann et al. 2002). These studies suggest that hybridization between introduced biological control agents needs further study as it may present an opportunity for true shifts in host use patterns (Szűcs et al., 2019a).
Most research on hybridization of biocontrol agents relies on F1 and sometimes F2 crosses generated in the lab (e.g. Bitume et al., 2017;Szűcs et al., 2012). However, the consequences of hybridization for both fitness and host use may not be evident within that time frame. Furthermore, later generation crosses or crosses among more than two species can occur when hybridization is ongoing in natural environments. Selection and drift can also act on the new hybrid combinations, increasing the frequency of some allelic combinations over others (Szűcs et al., 2011). To characterize fitness and host use of natural hybrids of complex ancestry, high-resolution estimates of ancestry (which species, or mix of species, are represented in an individual's genome) are needed for individuals sampled from the field. These estimates are now possible with modern genomic tools (Andrews et al., 2016;Stahlke et al., 2022).
Here, we study the consequences of hybridization among four species of Diorhabda (Coleoptera: Chrysomelidae), which are biocontrol agents with similar host preferences introduced to control invasive Tamarix spp. (Caryophyllales: Tamaricaceae) in North America (Milbrath and DeLoach, 2006a). Hybrids between three Diorhabda species are common in the introduced range (Stahlke et al., 2022), and previous studies indicate that shifts in host preference are possible (Thomas et al., 2010) and hybridization may play a role in host choice (Bitume et al., 2017). Two non-target plant species in the introduced range are palatable to Diorhabda and it is important to understand the potential impacts of biological control on them and whether that risk has changed since original host testing. We ask the following questions: 1) In locations where species overlap, what hybrid individuals (with mixed ancestry) are found? 2) How does hybridization affect fitness proxies of body size and fecundity? 3) How does hybridization affect palatability and adult feeding preferences on target and non-target hosts? 4) How does hybridization affect larval performance on non-target hosts? We synthesize these results to assess the role of hybridization in risk and efficacy for Tamarix biocontrol.

Study system
The Tamarix-Diorhabda biocontrol system provides a unique opportunity to examine the effects of hybridization. Four species in the genus Diorhabda were released in the United States. Hybrids between D. carinata (Faldermann), D. elongata (Brullé), and D. sublineata (Lucas) are viable and have been found in the field since 2013, or at least 16 generations (Bean et al., 2013;Knutson et al., 2019;Stahlke et al., 2022). The fourth species, D. carinulata (Desbrochers), does not readily hybridize (Bean et al., 2013). The four species have very similar external morphology and are only distinguishable through dissection (Tracy and Robbins, 2009) or genetic analysis (Bean et al., 2013), the latter of which was done in this study.
All four Diorhabda species have similar host preferences (Milbrath and DeLoach, 2006a) and survival is highest and oviposition is preferred on Tamarix chinensis Lour. × Tamarix ramosissima Ledeb. hybrids (DeLoach et al., 2003;Lewis et al., 2003a;Milbrath and DeLoach, 2006a). T. chinensis × T. ramosissima hybrids (hereafter, target host or tamarisk) are the most common genotype of invasive tamarisk in North America (Gaskin and Kazmer, 2009), but introgression between the two species varies by latitude, with T. ramosissima dominant in northern areas and T. chinensis dominant in southern areas (Gaskin and Kazmer, 2009;Williams et al., 2014). Tolerance of and resistance to herbivory also varies with tamarisk introgression, with T. ramosissima being more tolerant of and less resistant to herbivory (Williams et al., 2014).

Insect and plant material
Diorhabda were sampled off of the target host (tamarisk) at ten sites across eastern New Mexico and western Texas, USA from 16 to 19 September 2019 (Fig. 1, Supplemental Table S1, sites labeled A-J from north to south), a region where hybridization between the three species had been confirmed (RiversEdge West, 2021;Stahlke et al., 2022). At seven sites, at least 20 adults were collected for direct genotyping. The remaining three sites were chosen for both phenotyping and genotyping, as previous research suggested they should contain hybrids across a gradient of ancestries (Stahlke et al., 2022). From these three sites, 200 adult Diorhabda were collected, brought to the lab, and maintained in growth chambers with a 16:8hr light:dark regime at 28 • C : 20 • C day: night (Bean et al., 2007). Diorhabda were allowed to reproduce on the target host for one generation in the lab to reduce maternal effects prior to collecting phenotypic data on correlates of fitness and host use. At the time of collection and for all experiments, ancestry, including hybrid status, was unknown for all individuals, as those data could only be collected through lethal sampling after experiments ended.
The target host plant (T. chinensis × T. ramosissima) was collected as vegetative cuttings from Bonny Reservoir, Colorado, USA in March 2019. Since the level of introgression between T. chinensis and T. ramosissima impacts Diorhabda herbivory (Williams et al., 2014), we genotyped many cuttings and used only cuttings that were 30-44 % T. ramosissima for all rearing and host use testing. This tamarisk genotype was chosen as it represents tamarisk ancestry commonly found in the hybrid zone. AFLP markers were used to determine the level of introgression of each cutting following Gaskin et al. (2012). Genomic DNA was extracted from approximately 20 mg of silica-dried material using a modified CTAB method (Hillis et al., 1996). The AFLP method followed Vos et al. (1995) with modifications as in Gaskin et al. (2012) using the two polymorphic primer pairs MseI + CAT/EcoRI + ACC and MseI + CTA/EcoRI + ACC. Loci were initially scored by the fragment analyzer software GeneMapper v. 4.0 (Applied Biosystems). These bins were then manually screened, making this a semi-automatic scoring method, as suggested by Papa et al. (2005). Average estimated admixture coefficient (Q, or assignment value) of each plant to T. chinensis and T. ramosissima was determined using Structure v. 2.3.3 (Falush et al., 2007(Falush et al., , 2003Pritchard et al., 2000). Prior population information was included for 100 native Asian Tamarix plants included for reference (but not the plants collected from the USA), admixture was assumed to be possible, allele frequency was considered to be independent in each population and 50,000 burn-in and 50,000 run lengths were used.
Athel was collected near Yuma, Arizona, USA and Frankenia was collected near Lone Pine, California and Goleta, California, USA and neither were standardized by genotype. All three host species were grown in the Insectary Greenhouse at Colorado State University during the experiments and all plants were watered and fertilized regularly.

Genomic analysis of ancestry
We generated genomic data to determine ancestry for 20 adult Diorhabda collected from 7 sites and all female adults from 3 additional sites that were phenotyped in choice and no-choice host use trials (see section 2.5). As noted above, ancestry was unknown at the time of the host use trials, because our genomic approach used whole head and thoraxes, and thus any phenotyping had to occur prior to DNA extraction. All individuals were frozen at − 80 • C either after field collection or after completion of host use trials to protect high molecular weight DNA. Genomic data and ancestry analyses were performed following Stahlke et al. (2022). Briefly, DNA was extracted from the heads and thoraxes of individuals using a Qiagen DNeasy Blood and Tissue Kit following the manufacturer's protocol and treated with 4 μL RNase A (Qiagen) to eliminate RNA contamination. DNA sample concentration was quantified for each individual by fluorometric quantification (Qubit 2.0 HS DNA assay; Invitrogen, Life Technologies). We prepared a total of 321 individually barcoded RADseq samples, including two individual replicates per plate, across four single-digest RADseq libraries using the 8-bp restriction enzyme SbfI following the protocol described by Ali et al. (2016). Adapter-ligated libraries were multiplexed to achieve approximately 47.83 million 150 bp paired-end reads per library sequenced on an Illumina NovaSeq 6000 (University of Oregon).
Genotyping, filtering, and Bayesian clustering analyses were generally performed as previously described (Cerca et al., 2021;Rochette et al., 2019;Stahlke et al., 2022). For these analyses we included 40 RADseq samples previously identified as pure individuals of each of the four Diorhabda species derived from previous work to serve as a reference panel (Ravagni et al., 2021;Stahlke et al., 2022). This reference panel included individuals from all source populations introduced to North America, sampled from both laboratory cultures and field-caught individuals from the native and introduced ranges (data available at: Clark et al., 2022b). Briefly, all samples were aligned to the D. carinulata reference genome (Genbank Accession GCA 020975425.1) with bowtie2 version 2.2.9 (Langmead and Salzberg, 2012) and sorted with SAMtools version 1.9 (Li et al., 2009). Genotyping was performed including the reference panel with gstacks, using the default 'maruki_low' model, which incorporates population allele frequencies (Maruki and Lynch, 2017). We performed an iterative procedure to mitigate the effect of allele dropout and missing data within and among populations (Cerca et al., 2021). Through this process, samples with < 4x effective coverage, as well as individuals or loci with>50 % missing data within and across populations were discarded from further analysis.
Using a single SNP from each RAD locus and the program Structure version 2.3.4 (Pritchard et al., 2000), we identified species-specific ancestry clusters based on the reference panel of known single-species individuals. As previously validated (Stahlke et al., 2022), we used the uncorrelated allele frequency model and allowed the alpha parameter to be inferred for each population (Falush et al., 2003). We executed 10 independent runs for each K from 1 to 10,allowing a burn-in period of 10,000 steps and 10,000 Markov chain Monte Carlo replicates, and  Table S1. Unfilled points represent sites from which first-lab generation adults were genotyped following phenotyping. Most of the larval performance data came from site A. Filled points represent sites from which field individuals were genotyped.
printed the estimation of 90 % credible intervals. After inspection of likelihood values for each K, we used the clustering results of K = 4, corresponding to each of the four species, to infer ancestry across all field-collected and experimental individuals (Stahlke et al., 2022) (Supplemental Fig. S1). We examined the confidence intervals across independent runs to conservatively identify the threshold at which operational pure ancestry could be confidently inferred for diagnostic individuals for each species, q = 0.024 (i.e., the lower bound of the 90 % credible interval), below which pure identification could be unreliable and due to technical biases (Caniglia et al., 2020). Thus, by including reference samples of known identity, individuals with ancestry from more than one species (mixed ancestry) could be identified as having a hybrid origin resulting from secondary contact in North America.

Fitness
Body size and fecundity were used as fitness proxies. For insects in general (Berger et al., 2012) and Diorhabda specifically (Clark et al., 2022a), body size is a good indicator of fitness since it is associated with competitive ability and fecundity. These fitness proxies were measured on individuals reared in a common environment in the lab and thus differences among individuals from different collection sites or with different ancestry indicate evolutionary shifts and not plastic effects caused by differing field conditions. Adults from the first lab generation were weighed at emergence. Eggs laid during the 24-hour choice test (see section 2.5) were counted as a measure of fecundity.

Palatability of host plants and preference of adults
We measured palatability of each host plant through no-choice tests in which adults were presented with only one host plant and feeding was measured. We measured adult feeding preference through choice tests between all three host plants, in which adults could choose which plant to feed when presented with multiple options. All females were paired with a male from the same collection site shortly after eclosion and kept in pairs until the host use trials. After all females were laying eggs (about four days after emergence), each female was randomly assigned to a nochoice test with a single host plant or a choice test between all three host plants. Trials were started daily for five days from 13 to 17 November 2021. Choice and no-choice tests were similar to those conducted by DeLoach et al. (2003) and are briefly described here. Bouquets of each plant were made from 10 to 15 cm cuttings of fresh foliage with the cut end in a vial of water. For no-choice tests, a single bouquet of one plant was placed on a 5x10 cm piece of weighing paper at the center of a 15.25 cm diameter petri dish. In choice tests, separate pieces of weighing paper were evenly spaced in a petri dish and a bouquet of each plant was placed on each paper. In both types of tests, a single female beetle was placed at the center of the dish, away from any one plant, and allowed to eat and oviposit for 24 h. We quantified feeding behavior by measuring frass (insect excrement) left on and under each plant after 24 h (Bitume et al., 2017;DeLoach et al., 2003). At the end of 24 h, frass was swept off each bouquet with a paintbrush onto the corresponding weighing paper and any plant material was removed from the paper. For each trial, the weighing paper with frass was photographed with a ruler for scale. Photographs were analyzed by first reducing shadows in Windows Live Photo Gallery editor to improve photo quality, then area of frass under each plant (in mm 2 ) was determined with the Analyze Particles function in ImageJ (Abràmoff et al., 2004). All females from no-choice and choice trials were frozen at − 80 • C for genomic analysis after completion of the trials.

Performance of larvae
Eggs were collected from the adults that were assigned to the choice test prior to the 24-hour trial and during the choice test if eggs were laid on the target host. Due to space and foliage constraints, most replicates came from site A (n = 32), with only a subset of individuals from sites F (n = 6) and I (n = 4). On the day larvae hatched, larvae from each fullsibling family were randomly assigned in groups of 6-10 to develop on a single host plant (target host, athel, or Frankenia). Larvae from each fullsibling family were assigned to all three hosts as they hatched, over a 7day period from 20 to 26 November 2021. Bouquets of each plant were replaced every 1-2 days, as needed. Twelve days after hatching, the number of surviving larvae in each container was recorded. All surviving larvae were frozen at − 20 • C and then dried at 70 • C for 24 h before weighing each individual to the nearest 0.01 mg. Survival and weight were recorded at 12 days in order to obtain a standardized measure of development before the fastest growing larvae reached the pupa stage.
During the experiments, an infestation of leaf hoppers and spider mites attacked the target host in the greenhouse, so availability of highquality target host material decreased during the larval performance tests. Since plant quality can influence host performance (Lewis et al., 2003a), the data from larval performance on the target host have been excluded from analyses since they reflect low plant quality, not true performance. As our scientific question focuses on performance on nontarget hosts, this exclusion restricts our ability to compare to the target host, but does not compromise our ability to address our question. Larvae were not genotyped.

Statistical analysis
To test for relationships between ancestry and individual phenotypes of fitness and host use, ancestry was described with categorical scores and continuously. The categorical scoring of genotypes among the phenotyped individuals produced six categories of hybridization: pure D. carinulata (where "pure" here indicates q > 0.976 for a single species, based on the confidence of the ancestry assignment), pure D. carinata, pure D. sublineata, hybrids of D. carinata × D. sublineata (where "hybrid" here indicates q > 0.024 for more than one species), hybrids of D. carinata × D. elongata, and hybrids of D. carinata × D. sublineata × D. elongata. We used this broad classification of pure species and hybrids that captures any amount of introgression above the 0.024 threshold because we are interested in the current field populations rather than in early generation crosses (which are reported on in Bitume et al. 2017). To ensure the classification threshold did not impact the results of the host specificity analyses, we classified hybrids using a more conservative q > 0.1 threshold, and found no qualitative differences in the results of any analysis. Additionally, we visualized individual genotypes with principal components analysis (PCA) using adegenet version 2.1.8 (Jombart, 2008) as a complementary inspection of hybridization classifications. We found that individuals identified as hybrids clustered together, between the two parental species, indicating that our classification of hybrid categories is robust (Supplemental Fig. S2B).
Hybrids of D. carinata × D. elongata were excluded from further analyses due to very low realized sample sizes (n = 4). Pure D. carinulata beetles were also excluded because they rarely hybridize with the other species and are not relevant to the questions of this study. Thus, following genotyping, we could compare phenotyping results from two pure species (D. carinata and D. sublineata), D. carinata × D. sublineata hybrids (hereafter, 2-way hybrids), and D. carinata × D. sublineata × D. elongata hybrids (hereafter, 3-way hybrids).
Ancestry was also described continuously. Proportion D. sublineata and D. carinata ancestry were tightly correlated because most individuals in our samples were either D. sublineata, D. carinata, or a hybrid of the two. Proportion D. carinata ancestry was chosen for the statistical models because it was a slightly better predictor in most models and in no models were results qualitatively different between the two.
Adult weight was analyzed using models with Gaussian error distribution, with ancestry (separate models for categorical (factor with 4 levels) and continuous (proportion D. carinata) ancestry) and collection site as fixed effects, and emergence date as a random effect. Both ancestry and collection site were included in the model in order to estimate the effect of each after accounting for the effect of the other. For fecundity, models with negative binomial error distribution and zeroinflation were fit to egg counts, with ancestry (categorical and continuous) and collection site as fixed effects, weight as a fixed covariate, and test date as a random effect. The significance of model effects here and below were tested with Wald χ 2 tests with degrees of freedom denoted by subscripts. Pairwise comparison of means was performed with t-tests with degrees of freedom denoted by subscripts in the emmeans package (Lenth, 2020).
Palatability of each plant was estimated from no-choice tests. A mixed model was fit to frass area (mm 2 ) with gaussian error distribution with ancestry (categorical and continuous), plant treatment (3 levels: target host, athel, Frankenia), and their interactions as main effects, weight as a covariate, and trial date as a random effect. Sample size per ancestry group per treatment ranged from 3 to 13 (average = 7). Additional models with collection site as a fixed or random effect were fit and compared using AIC.
Feeding preference for each host plant was estimated from the choice tests. We modeled total frass area under each plant and proportion of total frass under each plant separately. For frass area, a small constant of 0.01 was added to the frass area under each plant in a choice arena before log transformation in order to fit the assumptions of a linear model. Models with gaussian error distribution were fit with ancestry (categorical and continuous), plant in the choice arena, and their interaction as main effects, weight as a covariate, and date of trial and individual beetle as random effects. We also fit additional models with site as a fixed or random predictor. To account for differences in variance across predictors, additional models were fit that included a dispersion component (Brooks et al., 2017). AIC was used to compare between models and results are reported from the model with the lowest AIC. The proportion of frass under each plant was calculated by dividing frass area under one plant by the total frass area for each individual. Proportion frass was then transformed with the formula, (x(n − 1) + 0.5)/n, where x is the proportion and n is the sample size, so the data were between 0 and 1 to fit the assumptions of a beta regression (Douma and Weedon, 2019). A separate beta regression model was fit for each of the three plants, with ancestry (categorical and continuous) as a main effect, weight as a covariate, and trial date as a random effect. Additional models with dispersion components and collection site effects were fit and compared using AIC.
For larval survival, a model with binomial error distribution was fit to proportion survival for each family, with ancestry of mother (categorical and continuous), plant (two levels: athel, Frankenia), and all interactions as main effects, weight, date of hatching (correlated with plant quality), and collection site as fixed covariates, and mother as a random variable, to account for the nested design. For larval growth, a model with gaussian error distribution was fit to dry weight of larvae at 12 days of development with ancestry of mother (categorical and continuous), plant, and all interactions as main effects, and weight and hatch date as fixed covariates. To account for the design of the experiment, mother and the interaction between mother and treatment were random effects. Additional models with dispersion components were fit and compared with AIC.
Host specificity can be reinforced by adult individuals preferring hosts that increase the performance of future generations (Poisot et al., 2011). This preference-performance association can be an important indicator for how host specificity may evolve. We tested this association by calcuating Pearson correlation statistics between mother's preference (proportion of total frass under one plant in a choice test) and larval performance (proportion survival on one plant).

Ancestry of Diorhabda
We determined the ancestry composition of at least 20 individuals per site from the field from seven sites across the region and 59 to 63 per site females from the first lab generation from three additional sites that were also phenotyped in the lab. After read-processing and aligning to the existing reference genome for D. carinulata (Stahlke et al., 2022), effective per-sample coverage across the 38,641 loci had a mean of 40x and standard deviation of 20.4x. Structure analysis was performed on 1,838 SNPs. At K = 4, change in likelihood values plateaued (i.e., the Ln''(ΔK) method; Supplemental Fig. S1) and matched species identities for the reference panel of four species. Structure results for K = 5-10 were the same as those for K = 4 (I.e., no additional clusters were recovered with increasing values of K), so we present only K = 4.We found three of the four possible parental species (finding no D. elongata) across our collection sites. Across the entire region, 25 % of individuals were assigned to be pure D. sublineata, 9 % pure D. carinata, and 3 % pure D. carinulata (D. carinulata is expected to be rare in the region sampled). 62 % of all individuals were assigned to be hybrids, with hybrids of D. carinata × D. sublineata dominating (69 % of hybrids) (Fig. 2). We found evidence of one hybrid with D. carinulata (out of 400 individuals genotyped), which is a cross that is expected to be rare (Bean et al., 2013;Stahlke et al., 2022). Ancestry composition differed between the three phenotyped sites (A, F, and I) ( Table 1, Supplemental Fig. S2A), but it was overall similar to the distribution of ancestry throughout the region (Fig. 2). For example, phenotyped site A is similar to sites B and D with D. carinata ancestry highly represented among the individuals. Phenotyped site F is similar to sites C, E, and G with D. sublineata highly represented. Phenotyped site I is similar to site H, with a high proportion of 3-way hybrids.

Fitness
Results of fitness and host-testing experiments are reported only for the four ancestry groups that had sufficient sample size to analyze: D. carinata, D. sublineata, 2-way hybrids (D. carinata × D. sublineata), and 3-way hybrids (D. carinata × D. sublineata × D. elongata). Average adult body size ranged from 9 to 22.5 mg. Collection site was a statistically significant predictor of size (χ 2 2 = 40.59, P < 0.001), while ancestry (categorical) was not (χ 2 3 = 2.96, P = 0.398). On average, when measured in a common lab environment, individuals from sites I (mean = 16.8 mg, SE = 0.38) and A (mean = 16 mg, SE = 0.37) were larger than individuals from site F (mean = 13.7 mg, SE = 0.364; pairwise tests P < 0.001), but not different from each other (pairwise test P > 0.05) (Fig. 3A). Within hybrid beetles, proportion D. carinata ancestry was not significantly associated with adult weight (χ 2 1 = 0.33, P = 0.566). Average fecundity for each ancestry group ranged from 17.8 to 20.0 eggs in 24 h. There were no significant differences between ancestry groups or collection sites in fecundity (Fig. 3B).

Palatability of host plants and preference of adults
No-choice tests allow us to assess the palatability of each host plant, or the amount an adult will feed when presented with no other options. In no-choice tests, the target host and athel received more frass than Frankenia (both pairwise P-values < 0.01), but there was no difference between the target host and athel (t 60 = 0.58, P = 0.833) (Fig. 4). Ancestry (categorical) nor its interaction with host plant were significantly associated with palatability (ancestry: χ 2 3 = 2.75, P = 0.432; interaction: χ 2 6 = 3.63, P = 0.726). Larger beetles left more frass than smaller beetles (χ 2 1 = 5.63, P = 0.018). Results were qualitatively the same when collection site of origin was included in the model.
Total frass area under each plant in a choice test allows us to assess preference and potential feeding damage for each host plant in absolute terms and to compare preference for host plants within each ancestry group. The best fitting model modeled dispersion of errors independently for each plant (see section 2.7 Statistical Analysis) and conclusions did not differ when collection site of origin was included in the model, so it was excluded. In general, Diorhabda preferred the target host over non-target athel and Frankenia, but the effect differed by ancestry (Fig. 5). D. sublineata, 3-way hybrids, and 2-way hybrids, all preferred the target host over both athel and Frankenia (all pairwise contrasts between plants P < 0.01). For D. carinata, there was no difference in preference between the target host and athel (t 231 = 0.48, P = 0.629), but both were preferred over Frankenia (P-values < 0.001).
Proportion frass under each plant for each trial in a choice test allows us to assess preference for each host plant relative to the other plants in the choice arena and estimate the effects of ancestry on preference for one host at a time. Including collection site of origin in any model was Fig. 2. Diorhabda ancestry assignment for individuals from 10 sites (A-J, from north to south) and samples of pure species used for reference. Each vertical bar represents a single individual and the proportion of their genome that was confidently assigned (q > 0.024) to each parental species. Sites are labeled as on Fig. 1 and Tables 1 and S1. Individuals from bolded sites (A, F, and I) were from the first lab generation and phenotyped in the lab. Individuals from the other 7 sites were collected from the field.

Table 1
Ancestry of Diorhabda in each phenotyped collection site. D. carinata × D. elongata hybrids and D. carinulata were excluded from fitness and host-use analyses due to low sample sizes and lack of hybridization, respectively. not preferred by AIC and did not change any conclusions, so it was excluded. Preference for the target host was above or not significantly different from 50 % for all ancestry groups, but D. carinata tended to have lower preference for the target host than D. sublineata (pairwise comparison t 76 = 1.98, P = 0.051) (Fig. 6A). D. carinata had preference for athel of 56 %, which was significantly stronger than D. sublineata's preference of 29 % for athel (pairwise comparison t 76 = 2.14, P = 0.035) (Fig. 6B)  way hybrids) had preference for Frankenia over 75 %, while no D. carinata or D. sublineata individuals had preference for Frankenia over 12 % (Fig. 6C). Higher proportion D. carinata ancestry was associated with a decrease in preference for the target host (χ 2 1 = 4.35, P = 0.037) (Fig. 6D). However, when pure species were excluded from the analysis, the trend was in the same direction, but no longer significant (χ 2 1 = 0.69, P = 0.405). This indicates that the trend is driven by the differences between pure D. sublineata and D. carinata, rather than introgression within the hybrids. Preference for athel tended to increase with proportion D. carinata ancestry, but the trend was not significant (χ 2 1 = 2.62, P = 0.106) (Fig. 6E). Preference for Frankenia was best modeled by accounting for uneven dispersion (ΔAIC = 38.14). Preference for Frankenia significantly increased with proportion D. carinata ancestry (χ 2 1 = 68.13, P < 0.001) (Fig. 6F). When pure species were excluded from this analysis, the trend was still highly significant (χ 2 1 = 14.61, P < 0.001), indicating that D. carinata ancestry within hybrids is associated with increased preference for Frankenia.

Performance of larvae
Performance to 12 days of development on the non-target hosts allows us to estimate the ability of larvae to feed and develop on the nontarget hosts, if adults oviposited on the non-targets. Larval survival was higher on athel than on Frankenia across all ancestry groups (athel mean = 0.929, SE = 0.0259; Frankenia mean = 0.566, SE = 0.059; χ 2 1 = 7.40, P = 0.007). Larval survival differed by ancestry of the mother on both non-target host plants. Larvae from 3-way hybrid mothers had the highest survival overall and it was significantly higher than survival of larvae from 2-way hybrid mothers (P < 0.05), marginally higher than survival of D. carinata mothers (P < 0.1), but not significantly different from survival of larvae from D. sublineata mothers (P = 0.26) on both athel and Frankenia (Fig. 7A). However, 3-way hybrids were represented by only three families in larval performance tests.
Weight of larvae after 12 days of development on a non-target plant was best modeled with additional dispersion parameters for ancestry of the mother, treatment, and their interaction. There was a significant interaction between ancestry of the mother and treatment (χ 2 3 = 15.06, P = 0.002), such that larval weight significantly differed by ancestry group on athel, but not on Frankenia (Fig. 7B). On athel, larvae from D. carinata and 2-way hybrid mothers were larger than larvae from D. sublineata and 3-way hybrid mothers. Weight of mothers (χ 2 1 = 1.52, P = 0.218) and hatch date of larvae (χ 2 1 = 0.01, P = 0.913) were not significantly associated with larval weight.

Preference-performance correlation
There was no significant trend between maternal preference and larval performance on either athel or Frankenia (Supplemental Fig. S3).

Discussion
Hybridization of biocontrol agents after introduction may have consequences for fitness and host use of hybrid individuals (Ellstrand and Schierenbeck, 2000;Szűcs et al., 2019b). Here, we examine the effects of multiple generations of hybridization on fecundity, body size, and host use of Diorhabda using both genomic and phenotypic data. Overall, we find that hybridization has largely not impacted body size or fecundity, and that host use of hybrid individuals is similar to that of the parents, with a few exceptions. We discuss our findings in the context of evolutionary theory and the implications for biocontrol of Tamarix and for biocontrol programs generally.
The host use of specialized phytophagous insects is often very stable (Hardy et al., 2020). Our results indicate that both pure species and hybrid Diorhabda prefer the target host, tamarisk, in choice tests of adult feeding preference. Individuals of pure D. sublineata ancestry chose the target host 68 % of the time, which matches well with previous studies, where preference for the target host ranged from 55 to 80 % with an average of 68 % across all studies (standardized for comparison between the target host and athel only) (Bitume et al., 2017;DeLoach, 2006a, 2006b) (previous studies summarized in Supplemental  Table S2 and Supplemental Fig. S4). In this study, preference for the target host of individuals of pure D. carinata ancestry was not significantly different from 50 %, somewhat lower than previous studies, which ranged from 54 to 71 %, with an average preference of 63 % (Bitume et al., 2017;DeLoach, 2006b, 2006a). The low preference for the target host by D. carinata resulted in a negative trend between D. carinata ancestry and preference for the target host and a positive trend for preference for athel. The differences in preference of D. carinata between this and previous studies could be explained by two non-mutually exclusive mechanisms. First, D. carinata could have evolved to be less host specific since previous testing. Second, quality or genotype of the target host used in this study could impact host preference. A spider mite outbreak on the target host during the choice tests could have altered concentrations of secondary compounds that D. carinata were more sensitive to than D. sublineata or hybrids. Plant genotype may also impact our results. We standardized the genotype of the target host plant to reflect 30-44 % T. ramosissima ancestry, which is representative of common genotypes in the hybrid zone, but is not as palatable as genotypes with more T. chinensis ancestry (Williams et al., 2014). Here, in this study, both types of beetle hybrids we tested chose the target host at a rate intermediate between both parental species, indicating that novel genetic combinations present in hybrids have not significantly altered preference for the target host. Our results are in agreement with previous host-testing studies showing that non-target effects are possible on both athel and Frankenia, and that hybridization does not worsen potential effects. Previous work has shown that hybrids of D. sublineata × D. elongata had reduced preference for tamarisk compared to pure parental species (Bitume et al., 2017). This cross was not detected in any of the ten sites surveyed in 2019, so risk of non-target effects from this hybrid combination is very low in this region.
Frankenia was the least preferred host plant by all ancestry groups. The parental species D. carinata and D. sublineata chose Frankenia on average 1-2 % of the time, which is very similar to original host testing of all Diorhabda species (Herr et al., 2014;Lewis et al., 2003a;Milbrath and DeLoach, 2006a;Moran et al., 2009) (Supplemental Table S2 and Supplemental Fig. S4). However, three individual adults deposited over 85 % of their frass on Frankenia, whereas all other individuals depos-ited<12 % of frass on Frankenia. Those three individuals were categorized as hybrids. One was assigned 71 % D. carinata, 7 % D. sublineata, and 22 % D. elongata, making it a three-way hybrid. The other two individuals were of predominately D. carinata ancestry (95-97.5 %), with some D. sublineata ancestry (2.5-5 %). For these three individuals, frass was observed under at least two plants, indicating they explored the arena during the 24-hour trial. Feeding on Frankenia increased with D. carinata ancestry in hybrid individuals. This indicates that hybrids may have genetic variation associated with feeding on Frankenia, and this variation could be selected on in the future. However, offspring of the individuals that preferred Frankenia performed better on athel than Frankenia and exhibited no increase in survival or growth on Frankenia relative to the offspring of mothers that did not prefer Frankenia. The lack of a correlation between maternal preference and larval performance suggests that preference for Frankenia will not be reinforced  (Gripenberg et al., 2010). Frankenia is primarily distributed in California, which does not currently overlap with the hybrid zone of Diorhabda. Based on the total evidence from this study, we believe that the risk to Frankenia from Diorhabda remains relatively low, but should be monitored in light of ongoing hybridization.
We found that all three host plants tested here were palatable in nochoice tests to adult Diorhabda, which was expected, based on previous host testing (DeLoach et al., 2003;Lewis et al., 2003a). Our results support previous work showing that feeding stimuli (e.g., plant volatiles) from Frankenia are reduced compared to the target host and athel (DeLoach et al., 2003;DeLoach, 2006b, 2006a). The utility of no-choice tests has been debated in the literature, since no-choice environments are rare in the field and may overestimate non-target impacts (Milbrath and DeLoach, 2006a;Schaffner, 2001). However, others have argued that these tests are useful, particularly in the Diorhabda system, because Diorhabda can quickly defoliate areas of the target host, leaving non-target hosts in what is essentially a no-choice situation in the field (Herr et al., 2009). Previous work on non-target impacts on Frankenia from Diorhabda show that in field situations, the risk to Frankenia is much lower than in artificial lab conditions, as were used in this study (Lewis et al., 2003a;Moran et al., 2009). Larval performance tests here only measured larval development to 12 days after hatching, which could over-estimate the probability of survival to adulthood. Our results show that risk to Frankenia already present from the biocontrol program is not increased by hybridization of Diorhabda, with a few exceptions noted above.
In a meta-analysis, hybridization among introduced species can increase fecundity and size, though these effects vary by taxon (Hovick and Whitney, 2014). Our results show that body size (which is related to fecundity and fitness in insects) measured in a common environment did not vary by ancestry. This is in partial agreement with a previous study showing no differences in fecundity or hatching rate between most labcreated hybrids and pure species of Diorhabda (Bitume et al., 2017). In that study, however, lab-created crosses between D. carinata and D. sublineata did exhibit some hybrid vigor (Bitume et al., 2017). The results of the present study do not support either fixed hybrid vigor or outbreeding depression in hybrids from the field, based on body mass and fecundity. Fecundity in this study was measured only over 24 h (during the choice tests), which likely explains the high proportion of females that did not lay any eggs. Diorhabda often lay eggs in clusters of 10-30 every 1-3 days (Lewis et al., 2003b). Body size did significantly vary by collection site when measured in a common lab environment, which indicates that there are genetically-based differences between collection sites. Collection site reflects multiple processes that may be occurring at each site, including introduction history of different species, population bottlenecks, and selection to different environmental factors. Future research will be needed to investigate this trend and determine if conditions present at each site or other factors are driving this pattern.
Secondary contact between species may have several long-term evolutionary outcomes, from complete isolation of species, to formation of stable hybrid zones, to complete admixture between them, depending on the isolation history of the hybridizing species, extent of genetic incompatibilities, and selection on the hybrid phenotypes in the local environment (Fischer et al., 2015;Sánchez-Guillén et al., 2016). While other studies have examined host use of lab-created crosses of biocontrol agents (Bitume et al., 2017), and change in proportions of hybrid biocontrol agents through time (Fischer et al., 2015), this study is the first to examine fitness and host use of a hybridizing biocontrol agent with hybrids collected from the field after many generations of hybridization. Our field-collected samples included beetles with pure ancestry of two species and variation in the amount of introgression between all three hybridizing Diorhabda species. These crosses are difficult to create in the lab but crucial for our understanding of adaptation across the range and the risk of host shifts in the field.
Hybrid and pure species beetles were found across the sampled region of New Mexico and Texas. Hybrids have been known to be in these locations since at least 2010 (Knutson et al., 2019), and the distribution of ancestry groups at each site have varied through time, indicating ongoing evolution across the region (Knutson et al., 2019;Stahlke et al., 2022). For example, in three sites that were sampled in this study (F, G, and H) and also by Stahlke et al. (2022), the abundance of individuals of pure ancestry has declined and abundance of hybrids has increased. At site F (Roswell, New Mexico), D. sublineata dominated from 2014 to 2017 (Stahlke et al., 2022), but in 2019 the site was primarily composed of D. carinata × D. sublineata hybrids. Similarly, at site G (Post, Texas), D. carinata was most abundant in 2014 (Stahlke et al., 2022), but in 2019, D. sublineata and D. carinata × D. sublineata hybrids were most abundant. Finally, at site H (Lake J. B. Thomas, Texas), D. elongata was previously abundant (Stahlke et al., 2022), but no D. elongata were found in 2019, and instead there were many 2-and 3way hybrids. Interestingly, although many releases of D. elongata have been made in the region (Deloach et al., 2012) and their establishment was confirmed both through genomic (Stahlke et al., 2022) and morphological (Knutson et al., 2019) characterization, we found no pure D. elongata individuals in 2019. We also found D. carinulata to be less abundant in the region than in previous studies and rarely forming hybrids with the other three species (Knutson et al., 2019;Stahlke et al., 2022). Introduction history, migration, genetic drift, and selection have likely all influenced the current distribution of hybrids through the landscape. A time series analysis of local ancestry and hybridization in these sites could further characterize evolution across the region (Gompert et al., 2017). Releasing multiple ecotypes or species of biocontrol agents that may be adapted to different environments in the native range has been suggested for many years to increase efficacy of an agent or control of the target pest through transient heterosis or increasing adaptive genetic variation (but see Clarke and Walter, 1995), though worries remain about hybridization leading to agents host-switching to non-target species (Szűcs et al., 2019b;Van Driesche and Bellows, 1996). Introgression has been beneficial for several biocontrol agents (Wright and Bennett, 2018), including the biocontrol agent of tansy ragwort Longitarsus jacobaeae (Szűcs et al., 2012), Dactylopius tomentosus, the biocontrol agent of cacti in the genus Cylindropuntia (Mathenge et al., 2010), and the egg parasitoid wasp Trichogramma chilonis Ischii (Benvenuto et al., 2012). Hybridization was detrimental to biocontrol efficacy of Dactylopius opuntiae, the biocontrol agent of cacti in the genus Opuntia (Hoffmann et al., 2002), but did not lead to host switching to non-target species. The present study corroborates previous work that shows that hybridization of Diorhabda seems to have no negative effect on fitness and may in some cases increase fitness (Bitume et al., 2017). It remains unknown how much hybridization may have facilitated longterm establishment of Diorhabda in the United States, but the widespread presence of hybrids over several years and increase in abundance of hybrids at some sites suggests that there may be selection for hybrids. Some have suggested that recent population contractions of Diorhabda in the hybrid zone could be due to hybridization (Knutson et al., 2019). Our results showing increased frequency of hybridization over time and similar fitness of hybrids and pure species do not support hybrid breakdown in this region. The ecological conditions that select for hybrid genotypes, their interaction with release history in the area, and consequences for establishment of agents and control of tamarisk need to be further explored (Stahlke et al., 2022).
In this study, we opted for a two-generation experiment of host use after only one generation in the lab to standardize maternal environment to most directly understand how contemporary field populations behave. Future work should explore if strong selection for feeding on Frankenia imposed over multiple generations will be able to increase preference for or performance on Frankenia, given our finding that some individuals prefer the non-target plant. Our study suggests that variation may be present for feeding and development on Frankenia, but that selection pressure may be lacking in the field. A multi-generation selection experiment would be able to assess rapid evolutionary change in host use and future risk under different scenarios of plant and insect range shifts (Müller-Schärer et al., 2020).

Conclusions
Hybridization between closely related species is hypothesized to be an avenue for increased genetic diversity that can influence fitness and host range of herbivorous insects. We find that fitness and host preference of hybrid Diorhabda biocontrol agents from the field is quite stable. We demonstrate the utility of genomic methods for understanding the effects of hybridization in the field between more than two species and the importance of using field-collected individuals to more fully understand the risks posed by hybridizing biological control agents. Our work supports the prediction and current evidence that hybridization among closely related species rarely generates transgressive phenotypes or incompatibilities, which suggests that hybridization is safe and even beneficial when the fundamental host range does not differ among parental species.

Data Availability
Raw sequence reads for each demultiplexed individual in this study are available through NCBI BioProject PRJNA902703. Ancestry assignment and phenotypic data are available at Ag Data Commons (https:// doi.org/10.15482/USDA.ADC/1528155).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.