Effect of mass rearing on the genetic diversity of the predatory mite Amblyseius swirskii

Amblyseius swirskii Athias‐Henriot (Acari: Phytoseiidae) is a predatory mite used to control whiteflies and thrips in protected crops. This biocontrol agent, originating from the Eastern Mediterranean region, has been mass‐reared for commercial use since 2005 and is widely used in augmentative biocontrol programs. As a polyphagous predator, it has to cope with different biotic and abiotic factors. However, possible adaptation to mass rearing for production might be hindering its resilience and capacity for optimum performance in the field. In this study, we investigated the effect of long‐term mass rearing on the genetic diversity of A. swirskii. We identified six microsatellite loci from whole‐genome nanopore sequencing of A. swirskii and used these in a comparative analysis of the genetic diversity and differentiation in eight wild populations collected from Israel in 2017 and a commercially available population. Our results indicate that the commercial population is 2.5× less heterozygous than the wild A. swirskii. Furthermore, the commercial population has the highest genetic differentiation from all the natural populations, as indicated by higher pairwise Fst values. Overall, we show that commercially reared A. swirskii have reduced genetic variation compared to their wild counterparts, which may reduce their performance when released to control pests in an integrated pest management (IPM) context.


Introduction
Food security and safety should be the cornerstone of modern agriculture and the overreliance on chemicals to control pests is jeopardising this goal. The accumulation of residues in the environment and the food chain, their negative impact on human health, biodiversity, soil and water quality, and the reduced efficacy due to resistance evolution have put the spotlight on biological control as the environmentally sound alternative for pest control. In augmentative biological control, mass-reared natural enemies are released in large numbers to protect the crops from the negative impact of certain pests (Eilenberg et al., 2001). Today, this approach is key for the success of integrated pest management (IPM) programs for many fruit and vegetable crops, such as cereals, maize, cotton, sugarcane, soybean, grapes, and many greenhouse crops on more than 30 million ha worldwide (van Lenteren et al., 2018).
The most widely used predatory mite in augmentative biological control is the polyphagous predator Amblyseius swirskii Athias-Henriot (Acari: Phytoseiidae), followed by two other phytoseiids, Phytoseiulus persimilis Athias-Henriot, and Neoseiulus californicus (McGregor) (Knapp et al., 2018). Amblyseius swirskii was originally described in 1962 from almond trees [Prunus dulcis (Mill.) DA Webb] in Israel (Athias Henriot, 1962), where it naturally occurs on various other annual and perennial crops, such as citrus, grapes, vegetables, and cotton (Swirski & Amitai, 1997). When whitefly resistance to pesticides caused outbreaks in greenhouses at the beginning of this century (e.g., Stansly et al., 2004), attention was drawn to A. swirskii because of its fast reproduction and high performance against the key pest Bemisia tabaci (Gennadius) (Nomikou et al., 2001;Calvo et al., 2015). Experiments showed that A. swirskii was able to suppress B. tabaci populations and also provided good control of thrips Frankliniella occidentalis (Pergande) on cucumber plants (Nomikou et al., 2002;van Houten et al., 2005). Today, this species is used to control whiteflies and thrips in greenhouse vegetables, fruits, and ornamentals in various parts of the world (Calvo et al., 2015). The importance of A. swirskii in current agriculture is showcased in protected sweet pepper crops in southeastern Spain, where the successful integration of this predatory mite in the IPM strategy against whiteflies and thrips has led to a sharp decrease in the use of chemical pesticides (Calvo et al., 2011;van Lenteren et al., 2018). In addition to its high efficacy in managing whitefly and thrips infestations, the success of A. swirskii has been further encouraged by its early establishment on crops before the target pests arrive, using pollen or factitious prey as a food source, and by its mass rearing on a plantless system using the stored-product mite Carpoglyphus lactis (L.) as factitious prey (Bolckmans & van Houten, 2006;Calvo et al., 2015).
The effects of long-term mass rearing on the quality and performance of commercial biocontrol agents such as A. swirskii in natural or semi-natural environments are largely unknown (e.g., Tayeh et al., 2012;Guzm an-Larralde et al., 2014;Rasmussen et al., 2018). Selection, inbreeding, and random genetic drift may lead to loss of genetic variability, loss of fitness, and reduced field performance of the mass-reared biocontrol agents (Mackauer, 1976). The constant and regular commercial mass-rearing conditions, which are optimized to obtain the maximum number of individuals in the shortest possible time, differ considerably from the greenhouse biotic and abiotic conditions. Furthermore, under these artificial conditions, searching time for preying, mating, or dispersal is intentionally restricted to accelerate and synchronize the breeding batches as much as possible. Therefore, long-term laboratory rearing eventually leads to selection for the facility conditions (domestication), which in turn affects the reproductive or behavioural traits of the arthropods (Hoffmann & Ross, 2018). For example, the reproductive fitness of Drosophila melanogaster Meigen populations captive for 50 generations was reduced in small populations because of inbreeding depression and in large populations due to genetic adaptation (Woodworth et al., 2002). Although domestication is likely to play a key role in the rearing and success of biocontrol agents, only few studies on natural enemies have addressed this process.
For instance, fecundity and emergence of the egg parasitoid Trichogramma galloi (Zucchi) were reduced on the target pest, when the parasitoids were maintained on factitious compared to natural hosts (Bertin et al., 2017). Long-term laboratory rearing of the codling moth parasitoid, Mastrus ridens (Horstmann), resulted in lower genetic diversity and higher occurrence of diploid males, most likely as an effect of the reduced allelic diversity on a complementary sex determination locus (Retamal et al., 2016). In inbred lines of this parasitoid, the sex ratio of males to females and the proportion of diploid males was higher, and fewer daughters were produced compared to outbred lines (Zaviezo et al., 2018). Also, inbred populations of the mass-reared parasitoid Trichogramma pretiosum (Riley) produced fewer offspring than genetically variable populations (Guzm an-Larralde et al., 2014). A severe lack of genetic diversity and phenotypic differentiation in wing size and abdomen colour were found for one inbred laboratory colony of the crop pollinator Eristalis tenax (L.) compared with natural populations (Francuski et al., 2014). In another case, reproduction parameters and virulence of the entomopathogenic nematode Steinernema glaseri (Steiner) were reduced after laboratory adaptation took place (Stuart & Gaugler, 1996).
The performance of biocontrol agents under variable field conditions depends on its genetic variability, which is defined by the initial size, origin, and degree of inbreeding of the founder colony of a biocontrol agent (Mackauer, 1976). If the commercial population is established from a limited number of individuals, then the genetic variation is expected to be low and further reduced due to inbreeding. The genetic variation in commercial populations may be reduced by random genetic drift during rearing as well, such as splitting the colony for commercial distribution, or the occurrence of bottlenecks (Nunney, 2003). One example is the aphid parasitoid Diaeretiella rapae (M'Intosh), which experienced a significant founder effect and a strong reduction in the genetic diversity upon its introduction in Australia, compared to the populations from the Old World (Baker et al., 2003). Finally, six biofactory populations of the cereal stemborer parasitoid Cotesia flavipes (Cameron) in Brazil showed genetic differentiation among each other, resulting in new guidelines for the mixture of populations in the mass-rearing strategy to enhance the genetic variability (Freitas et al., 2018).
The commercial A. swirskii was originally collected in Israel and is being reared on a large scale without its natural prey since 2005 (YM van Houten, Koppert Biological Systems, The Netherlands, pers. comm.). The aim of this study was to investigate the population genetic structure of field A. swirskii from Israel, as well as the genetic variation in a commercial, long-term reared population of this biocontrol agent. For this, we developed microsatellite markers for A. swirskii, using a cost-effective, in-house sequencing process, the MinION Nanopore sequencer for whole-genome sequencing, which was applied and evaluated for the first time in sequencing predatory mites. Also, in our study, we explore the application of pooled microsatellite analysis to predatory mites, thereby offering a cost-effective solution to low individual DNA yields in these minute organisms (Skalski et al., 2006).

Populations and DNA preparation
For the whole-genome sequencing, an inbred population from the commercially available Swirski-Mite â (Koppert Biological Systems, Berkel en Rodenrijs, The Netherlands) was established. Five females and five males from the general population were allowed to reproduce for approximately 10 days in a rearing system with cattail pollen (Typha latifolia L.) as food source, where eggs were collected and isolated on a new rearing setup. Subsequently, five F1 females and five F1 males were used to start a new generation. This was repeated for a total of 10 generations. Over 200 individuals from the 10th generation of the inbred population were collected in 96% ethanol. DNA extraction was performed using the Qiagen MagAttract HMW DNA Kit (Qiagen, Hilden, Germany) according to the manual, whereas the final DNA product was eluted with nuclease-free water.
For microsatellite analysis, in total eight populations were sampled, from five locations across Israel, from various host crops (citrus, kiwi, sweet pepper, Jerusalem artichoke, cotton) ( Figure 1, Table 1). Mites were collected from the field and transferred to the laboratory where isolated colonies were established and named based on the plant and location found, with Carpobrotus edulis (L.) NE Brown pollen as food source. After approximately 18 generations (6 months), 100 female mites were collected from each colony and placed in 96% ethanol at À20°C for DNA extraction. A single batch of 100 A. swirskii individuals from a commercially available Swirski-Mite bottle (Koppert) was collected directly and placed in 96% ethanol at À20°C. The commercial population was not reared in the laboratory before sampling, because this could cause a reduction in its genetic variation. Due to the small size of the mites (0.5 mm) and to reduce genotyping costs (Skalski et al., 2006), we performed pooled DNA extractions for individuals from the same population. Tissue homogenization was achieved by placing dry snap-frozen samples with beads in a shaker for 20 s. Finally, DNA was extracted using a high-salt extraction protocol (Maniatis et al., 1982).

Genomic sequencing, markers, and experimental procedures
The inbred A. swirskii population was sequenced using Oxford Nanopore Technologies (ONT) MinION sequencer (Ligation Sequencing Kit v.108, SpotON Flow Cell FLO-MIN107 R9.5) (ONT, Oxford, UK). Approximately 30 000 reads were generated through direct base calling with the MinKNOW platform (v.1.10.16) or offline base calling using Albacore (v.2.1.3) (ONT). All raw and uncorrected read files (both FAST5 and FASTQ formats) are available for download (Ferguson, 2018). The basecalled reads, in FASTQ format, were initially trimmed in CLC Genomics Workbench v.11.0 (Qiagen, Hilden, Germany) for a minimum read length of 100 bp with lowquality sequences removed (limit = 0.05), resulting in 11 511 trimmed reads (Qiagen). These trimmed FASTQ reads were then corrected using CANU v.1.4 pipeline (parameters: genomeSize = 100 m, correctedErrorRate = 0.120, stopOnReadQuality = false) (Koren et al., 2017). This resulted in a total of 574 corrected reads, all longer than 1 000 bp. These corrected, raw reads are deposited in the NCBI Sequence Read Archive (PRJNA433466). The corrected sequences were used for mining di-, tri-, tetra-, penta-, and hexanucleotide repeats using MSATCOMMANDER v.0.8.2 with the default settings, except for the number of repeats that was set to a minimum of three for tri-, tetra-, penta-, and hexanucleotides (Faircloth, 2008). The MSATCOMMANDER program uses Primer3 v.1.1.1 to design locus-specific, flanking primers. The settings used for Primer3 were as follows: no perfect repeats, product size 75-500, primer size: min 16, opt 20, max 24 bases, and primer annealing temperature: 56-64°C. Poor primers were removed based on the following criteria: duplicate, positive for hairpins, high complementarity, and selfing (Rozen & Skaletsky, 2000). The markers with their corresponding primer pairs were sorted by the count of repeats in descending order, as markers with fewer repeats can be less polymorphic and the first 24 primer pairs of microsatellite loci were selected to be tested for their efficiency using PCR.
A total of six microsatellite loci primer pairs amplified efficiently only the desired microsatellite marker and were selected for analysis, whereas the remaining were rejected for either low efficiency or low specificity. Primer specificity was considered low when the PCR resulted in products different from the expected size. The sequence data of the loci were deposited in the NCBI's GenBank with accession numbers MK267176-MK267181 (Table S1). The six primer pairs used for the microsatellite amplification and analysis were fluorescently labelled (Table S1). The microsatellite markers were amplified for the nine populations of A. swirskii according to the following protocol. The PCR reactions were performed in a total volume of 50 µl containing 60 ng of genomic DNA, 19 GoTaq polymerase buffer, 200 µM dNTP, 1.25 U GoTaq polymerase (Promega, Madison, WI, USA), and 0.4 µM of each primer. The PCR protocol was as follows: 5 min denaturation at 95°C, 35 cycles of [95°C for 30 s, 50-60°C for 30 s, 72°C for 60 s], followed by 10 min at 72°C. Primer sequences and annealing temperatures are indicated in Table S1. For analysis, microsatellite amplicons were diluted 2509 and pooled into sets of 2-3 in such a way that overlapping of alleles was avoided even if markers differed for the fluorescent dye (Table S1). The samples were denatured at 95°C and loaded with an internal size standard (GeneScan 500 LIZ) onto an ABI3730 capillary automated sequencer (Thermo Fisher Scientific, Waltham, MA, USA).

Microsatellite analysis
Genotyping of the pooled DNA samples was performed using the GeneMapper software (Thermo Fisher Scientific). Instead of allele frequencies, peak frequencies were calculated based on peak areas, which is more reliable than to correct for stutter bands in pooled microsatellite samples (Crooijmans et al., 1996). Microsatellite peak scoring was used to estimate allele frequencies as explained in Hillel et al. (2003). Peak frequencies lower than 0.05 were discarded and peak frequencies higher than 0.05 were subsequently recalculated to add up to 1 to correct for stutter bands that resulted from incomplete PCR cycles (Megens et al., 2008).

Data analysis
Neighbour-Joining clustering was performed for the microsatellite dataset, using Nei's D A genetic distance in POPTREEW with 1 000 bootstraps (Takezaki et al., 2014). Gene diversity was measured as heterozygosity for all loci and for all populations by using peak frequencies instead of allele frequencies according to Hillel et al. (2003). An average heterozygosity was calculated for each population across markers (H) and for each marker across populations (H m ). Similarly, the average allelic richness for each population across markers (A) and for each marker across populations (A m ), and the number of private alleles (N p ) for each population were calculated. Three genetic distance measures for various evolutionary models were calculated based on peak frequencies: Nei's genetic distance (D A ), Cavalli-Sforza chord measure (D CS ), and Reynolds' genetic distance (D R ) (Cavalli-Sforza & Edwards, 1967;Nei et al., 1983;Reynolds et al., 1983). Pairwise distances between each pair of the nine populations (81 estimates) were estimated for each measure. D A was chosen to be used in the downstream analysis as more appropriate, because this measure assumes that genetic differences are caused by mutations and genetic drift, contrary to the other measures, which do not take into account mutation (Cavalli-Sforza & Edwards, 1967;Nei et al., 1983;Reynolds , 1983). The distance matrices were used for multidimensional scaling (Gower, 1966) in R v.3.4.3 (R Core Team, 2018). Isolation-by-distance analysis with a Mantel (1967) correlation test based on 9 999 replicates was also performed using the 'ade4' package (Dray et al., 2007) in R v.3.4.3. Last, F st measures of genetic differentiation were calculated manually according to Wright (1984).

Results
Amblyseius swirskii genomic sequencing with the MinION generated 256 Mb total sequence data after base calling, or approximately 30 000 reads. After sequence trimming and correction, 574 corrected reads of at least 1 000 bp long where retained (5.2 Mb). Microsatellite mining of the genome sequences identified 2 423 microsatellite loci on 532 reads. The density of microsatellites in the genome is 466 microsatellites per Mb and 92.6% of our reads contain microsatellite loci. The most common type of microsatellites are dinucleotide repeats (72.1%), followed by trinucleotide repeats (23.7%), whereas the least common microsatellites found were the hexanucleotide repeats (0.12%) (Table S2). Primers designed using Primer3 were filtered for quality, resulting in 622 valid primer pairs. Of those, 24 primer pairs of the markers with the largest repeat count were selected and then tested for their specificity and efficiency. Finally, a set of six microsatellite markers was found to be specific and efficient for using in the microsatellite analysis of population pools. The nucleotide motifs of these six microsatellite markers are: di-(Asw1), tri-(Asw3, Asw6), and tetranucleotide repeats (Asw2, Asw4, Asw5) (Table S1). In total, 32 different alleles were scored across the six microsatellite loci for the nine populations, with four of those alleles being private, exclusively found in equal number of wild populations (Table 2). All markers were polymorphic in at least four of the nine populations and all of the populations had at least three markers polymorphic ( Table 2). The average gene diversity H m for the six markers analysed was 0.42 and the mean allelic richness per locus A m was 5.3 across populations and 2.9 within populations ( Table 3). The most polymorphic marker was Asw1 with 12 alleles across populations and a mean allelic richness of A m = 5.8 per population (Table 3). The gene diversity H m of marker Asw1 was 0.66 and all of the populations were polymorphic for this marker (Table 3). On the other hand, the least polymorphic marker was Asw5, with two alleles over all populations, with an average A m of 1.6 alleles per population and gene diversity H m = 0.19, and it was polymorphic in five of the nine populations (Table 3). No bias of field sample size on the allelic richness was found (Figure 2).
When comparing populations, the lowest frequency of polymorphic markers P, average heterozygosity H, and allelic richness A across loci, were found in the commercial population KBS-9 (P = 0.5, H = 0.21, A = 1.8). The average heterozygosity H in KBS-9 was 1.6-2.59 lower than the average heterozygosity present in wild A. swirskii (H = 0.34-0.53, P = 0.83-1, A = 2.5-3.7) ( Table 2, Figure 3). The heterozygosity H m of the commercial population KBS-9 ranged from 0.65 and A m = 4 alleles for marker Asw1, to 0 for markers Asw2, Asw3, and Asw5 (Table 3). The highest average heterozygosity across loci was found in A. swirskii population Had-4 (H = 0.53), collected from sweet pepper crops in Israel, and ranged from H m = 0.65 and A m = 3 alleles for marker Asw4, to H m = 0.38 and A m = 2 alleles for marker Asw2 (Table 3). The commercial population KBS-9 had the highest genetic differentiation from all the natural populations, as indicated by the pairwise F st values (Table 4).
The genetic distance measures Nei's D A , Reynolds' D R , and Cavalli-Sforza chord D CS were calculated between   populations and the highest average distances were found for the commercial population and for the population collected on kiwi (Table 2). Multidimensional scaling and Neighbour-joining clustering of the Nei's D A genetic distances derived from the microsatellite genotyping of the wild A. swirskii populations suggest three possible clusters, one is the population BHE-2-Kiwi separated from the other populations, a second including the populations Had-6-Citrus, HAs-7-Cotton, and HAs-8-Citrus, and last a cluster of the populations Had-5-Jerusalem artichoke, Had-4-Sweet pepper, HV-1-Citrus, and EY-3-Citrusthe commercial population KBS-9 falls in the latter cluster of wild populations (Figures 4 and 5).
The Neighbour-Joining analysis yielded results very similar to the Multidimensional scaling; however, moderate bootstrap values (<70%) were computed, because the number of microsatellite loci employed in the analysis was lower than the number of populations analysed ( Figure 5). No evidence for isolation-by-distance was found among the wild A. swirskii populations, according to the Mantel correlation test (simulated P = 0.30) and the linear regression fit (Figure 6).

Discussion
The objective of our study was to estimate and compare the genetic diversity in natural populations and a massreared commercial population of the biocontrol agent A. swirskii by analysing microsatellite markers. In the absence of a published genome of A. swirskii, the identification of suitable markers for the genetic population inference required the de novo sequencing of the genome. However, next-generation sequencing approaches, such as 'whole genome shotgun' (WGS) sequencing, need high DNA yield from a highly inbred population, which is intricate to obtain for this minute predatory mite because of its size. In this study, we selected the MinION sequencing platform because it allows small DNA input, and generates long sequence reads (Goodwin et al., 2015). Furthermore, the entire sequencing process, including library preparation can be done in-house and is cost-effective.
Despite the benefits of nanopore sequencing in generating genomic sequences from the inbred A. swirskii strain, the 30 000 sequences (256 Mb) were not enough to generate a genome assembly, because these sequences did not provide an adequate coverage of the whole genome. We therefore refrained from further assembling the reads into scaffolds, and used the raw reads as input for our microsatellite analysis. However, these reads were beneficial to mining and developing microsatellites for population genetic studies, and add to the low amount of publicly available sequence data for A. swirskii, and to the  microsatellite panels developed for phytoseiids of economic importance (e.g., Sabater-Muñoz et al., 2012). Microsatellite mining of the raw reads, primer design, and filtering yielded 2 423 microsatellite markers and 622 candidate primer pairs for population genetic analysis. However, when 24 of those primers were tested for efficiency and specificity, only a small fraction (25%) complied with both criteria. One explanation for the low success rate of our microsatellites could be the high error rate of the nanopore sequencing platform with 10% of the bases called wrongly (Laver et al., 2015) combined with insufficient read depth to properly correct for these errors. Another issue is the likely contamination of the sequencing starting material by the food source cattail pollen (T. latifolia). Contamination by the food source may have affected the outcome of the genome sequencing of A. swirskii providing sequence reads of the contaminant organism (in this case the T. latifolia) as well, and has been outlined as a putative hazard for predatory mite studies before (Hoy et al., 2003). Predatory mites are tiny organisms whose gut occupies most of their total volume and it is very difficult to be dissected and excluded, as it is the common technique in larger arthropod species studies. Also, the detection window of prey DNA in the guts of arthropod predators (analysed by PCR) can range from a few hours to various days post-feeding (King et al., 2008). Unfortunately, the genome of T. latifolia is not known yet and its sequences cannot be filtered out from the consensus genomic sequences. Hence, in order to improve the methodology, the food source provided to the predatory mites should be a species with a genome sequence available, as those sequences can be filtered out and removed from the sequencing reads. Furthermore, the coverage of  Individual genotyping of mites for population genetic inferences is an expensive and elaborate method, mainly because of their tiny size and the low amount of DNA available for extraction, which limits the number of loci analysed per individual. However, DNA pooling for assessing relative differences of allele frequency among populations can overcome these limitations and provide a costeffective alternative (Skalski et al., 2006). Although the effective number of peaks obtained from DNA pools is systematically overestimated compared to the actual number of alleles obtained from individual typing, those estimates can still detect relative differences in allele frequencies among DNA pools (Skalski et al., 2006). Microsatellitebased genetic diversity estimates obtained from DNA pools were found to strongly correlate with those obtained from individually typing in previous studies in chickens and pigs (Hillel et al., 2003;Megens et al., 2008). In those studies, the effective number of peaks correlated strongly with the effective number of alleles, as did the respective heterozygosity estimates. Hence, the use of DNA pools provides reliable estimates for the population's diversity (Hillel et al., 2003;Skalski et al., 2006;Megens et al., 2008). Also, the three genetic distance measures -D A , D CS , and D Rfor different evolutionary models correlate significantly between them, similar to that was found in chicken and pig lines for pooled samples, suggesting that the downstream inference based on D A is robust under alternative evolutionary scenarios as well (Hillel et al., 2003;Megens et al., 2008).
A small number of microsatellite loci (2 423) was identified in the genomic sequences of A. swirskii, similar to the low number of loci found in Tetranychus urticae Koch genome (Grbi c et al., 2011). However, the actual number of microsatellite loci is expected to be higher, because the genome of A. swirskii was partially covered, yet our estimate is useful for the calculation of the microsatellite density on the genome. For two other Acari speciesthe phytoseiid Amblyseius fallacis (Garman) and the ixodid tick Ixodes scapularis (Say)a low abundance of microsatellites has been proposed as well, following the generation of plasmid libraries enriched for microsatellites (Navajas et al., 1998;Fagerberg et al., 2001), although this method does not provide absolute estimates of microsatellite density to compare. In A. swirskii, dinucleotide repeats are the most abundant type of microsatellite, trinucleotides are less than half frequent, and longer repeats, such as tetra-, penta-, and hexanucleotides, are found markedly less, similar to most arthropod species (Pannebakker et al., 2010) and different from T. urticae, in which trinucleotides are the most abundant repeat motif (Grbi c et al., 2011). It is unclear why mites have a low microsatellite density compared to other arthropod species (e.g., Pannebakker et al., 2010;Abe & Pannebakker, 2017), but this does not need to reflect an overall low genomic diversity. For instance, despite the low microsatellite abundance, Van Zee et al. (2013) found a very high SNP density in the genome of the tick I. scapularis, suggesting SNPs to be a more suitable marker for population genetic analysis. Nevertheless, we did observe high levels of polymorphism in the six microsatellite markers employed in A. swirskii in the current study, as reflected by the finding that the least polymorphic marker was polymorphic in 50% of the wild populations and six out of eight markers were polymorphic in all of the wild populations. Our finding is in accordance with the elevated molecular evolution found for the phytoseiid mite Metaseiulus occidentalis (Nesbitt) (Hoy et al., 2016).
High genetic variation is observed in natural populations of A. swirskii, with the mean observed heterozygosity, H, ranging from 0.34 to 0.53, when compared to the commercial population (0.21), despite the fact that some natural populations were established by a low number of individuals. Heterozygosity estimates based on pooled microsatellite analysis found in our study are comparable to the average observed heterozygosity found in a field population of the predatory mite Neoseiulus womersleyi (Schicha) (Hinomoto et al., 2011). Considerable genetic variation was found also for tetranychid mites; T. urticae field populations from Europe had a spatial genetic structure along a latitudinal gradient (Carbonnelle et al., 2007). Tetranychus turkestani (Ugarov & Nikolskii) collected on crops and weeds in southern France had comparable average heterozygosity and mites living on different host plants did not demonstrate clear evidence for genetic differentiation, similar to our results (Bailly et al., 2004). Our results did not show a correlation between host plant and genetic diversity in the wild populations; however, the BHE-2 population collected on kiwi plants did show a lower genetic diversity, which might be due to the presence of high trichome density on this plant. Plant trichomes can affect negatively the predators by impeding their movement and causing their entrapment (Riddick & Simmons, 2014). Hence, the populations of predatory mites can be smaller on plants with high pubescence compared to glabrous plants. Natural populations are often characterized by higher levels of genetic variation compared to captive-bred organisms because genetic diversity is crucial for their survival under fluctuating environmental conditions and diverse ecosystems (Barrett & Schluter, 2008). For example, population genetic analysis of chicken using microsatellites showed that wild populations of the red junglefowl [Gallus gallus gallus (L.)] had 39 higher average heterozygote frequencies than a European, domesticated chicken population (Granevitze et al., 2007). Commercial A. swirskii were polymorphic only in 50% of the microsatellite markers and had lower average allelic richness and lower heterozygosity estimates compared to their wild counterparts, indicating a very low genetic diversity. Low genetic variation of a commercial biocontrol agent may affect their resilience when unexpected and even minor changes or fluctuations of the environmental conditions take place (Lommen et al., 2017;Wright & Bennett, 2018). For instance, in D. melanogaster it has been demonstrated that for higher inbreeding levels, the impact of environmental stress becomes significantly greater and can even lead to extinction (Bijlsma et al., 2000).
The decline in the genetic variation in the mass-reared population of this biocontrol agent may affect traits that are important in its performance, such as reproductive parameters, prey preference, or predation rate. Our laboratory colonies of wild A. swirskii populations were founded by low numbers of mites in some cases (8-100); however, still in all cases their genetic diversity was higher than in the commercial strain. Despite the very high numbers at which commercial A. swirskii populations are reared, genetic diversity can be low if the initial founding population is low (Mackauer, 1976;Bartlett, 1993). Furthermore, even if a large founding population was used, it may have included a few genotypes that are favoured under mass-rearing conditions, which can result in a bottleneck that reduces the genetic diversity (Nunney, 2003). As the field environment is different from the mass-rearing conditions, adaptation to the rearing environment can have a negative impact on the fitness of mass-reared biocontrol agents in the field (Sørensen et al., 2012). Genetic change in mass-reared populations has been described in many arthropods (Nunney, 2003), including experimental populations of D. melanogaster (Woodworth et al., 2002), but also in parasitoids reared for biological control (Guzm an-Larralde et al., 2014; Retamal et al., 2016;Bertin et al., 2017) and various insects reared for 'sterile insect technique' programmes since the 1980s (R€ ossler, 1975;Bush & Neck, 1976;Wong & Nakahara, 1978;Vargas & Carey, 1989).
To further improve the predatory mite A. swirskii for augmentative biological control and prevent the decrease in traits related to its fitness, we would suggest enhancing the low genetic diversity present in the mass-reared commercial A. swirskii by adding genetic variation from natural populations (Nunney, 2003). However, a better understanding of the levels of genetic variation in massreared A. swirskii populations should be obtained first, before new genetic material is added to mass rearing. The present study only tested a single commercial population limiting the scope of the conclusions. Additional genetic analysis of commercial populations can further assess the impact of genetic diversity on the performance of A. swirskii as a biocontrol agent, which remains one of the most successful natural enemies in augmentative biocontrol of whiteflies and thrips in greenhouse vegetables, fruits, and ornamentals in several parts of the world (Calvo et al., 2015;Knapp et al., 2018). Augmentative biological control of agricultural pests is the environmentally sustainable alternative to synthetic pesticides for safe food production and healthy ecosystems. We found that the genetic diversity of a widely used commercial biocontrol agent is limited compared to the natural conspecifics and this could have implications when the biocontrol agents need to adapt to adverse or novel conditions such as when being released to the field or greenhouse. Also, we showed that pooled microsatellite analysis provides a cost-effective method to determine the genetic diversity of minute biocontrol agents. Hence, adding to our knowledge approaches to determine the existing genetic variance of commercial biocontrol agents is imperative for the long-term success of IPM strategies and control programmes. Table S1. Primer pairs for the mitochondrial and microsatellite markers with corresponding four GenBank accession numbers, melting temperature (T m ), PCR product sizes, number of five repeats, motifs, and fluorescent dyes Table S2. The composition of the microsatellite repeats found in the genome of Amblyseius swirskii