DNA Barcoding and Phylogeny of Acari Species Based on ITS and COI Markers

C oxidase I) marker and ITS nuclear ribosomal DNA region for species identi ﬁ cation in Tetranychidae and Phytoseiidae. The molecular data allow us to establish species boundaries and phylogenetic relationships among several clades of Acari, mainly Tetranychidae and Phytoseiidae. Sequence comparisons between complete COI and the Acari mitochondrial COI, ITS1-5,8S-ITS2, and ITS2 among all Acari sequences have demonstrated that the selected regions, even small, gave enough informative positions for both species ’ identi ﬁ cation and phylogenetic studies. Analyses of both DNA regions have unveiled their use as species identi ﬁ cation characters, with special emphasis on Acari mitochondrial COI for Tetranychidae and Phytoseiidae species in comparison with the Folmer fragment, which has been universally used as a barcode marker. We demonstrated that the Acari mitochondrial COI region is also a suitable marker to establish a barcode dataset for Acari identi ﬁ cation. Our phylogenetic analyses are congruent with other recent works, showing that Acari is a monophyletic group, of which Astigmata, Ixodida, Mesostigmata, Oribatida, and Prostigmata are also monophyletic.


Introduction
Phytophagous and predatory mites, belonging to the Tetranychidae and Phytoseiidae families (Arachnida: Acari: Trombidiformes; Mesostigmata), respectively, are of economic importance in agriculture [1]. The first ones, tetranychids, because they damage the host plant, and even in some cases cause its death, meaning a reduction of production and economic value of crops [2]. The second ones, phytoseiids, because they render a positive value to growers and hence per establishment of multiple companies that are involved in their mass-rearing [3]. The control programs established against these plant pest species, which impact has been increased in the past years mainly due to the actual EU policy (2009/128/CE) for food safety and environmental protection, had required specific species identification [4,5], evaluation of economic thresholds for chemical treatments [6], and implementation of integrated pest management (IPM) procedures along the use of predatory phytoseiids as biological control agents [1,7].
As forehead mentioned, one key point in crop protection is the correct identification of the key pests and key natural enemies at the species level. In this sense, the taxonomy and systematics of Tetranychidae and Phytoseiidae have been traditionally based on morphology [8], but their minute size and reduced number of morphological taxonlinked structures make it a difficult task. These drawbacks have driven acarologists into the world of DNA markers in the past decades, which jointly with the use of sequence identity percentage and/or phylogenetic neighborhood allow species identification [4,[9][10][11]. Over the last years, several kinds of molecular markers have been used to study the Acari species with different purposes, from clarification of taxonomic problems [12][13][14][15][16][17][18], to the understanding of the dynamic of the pest [19], the establishment of dispersal patterns [19,20], the determination of population structure [21][22][23][24][25][26][27], their food preferences [9,11], and even to the determination of geographical origin of invasive species producing an outbreak [28][29][30][31][32]. Moreover, one species, Tetranychus urticae Koch (Trombidiformes: Tetranychidae), has been the focus of an International Genome consortium that has achieved the first full genome and transcriptome record for an economically important Acari agricultural pest [33]. Even having the whole genome of the model species available, and with the reduction of costs of Next Generation Sequencing (NGS), single-copy genes like mitochondrial Cytochrome C Oxidase I (coxI, here referred as COI) and the Internal Transcribed Spacer (ITS, here referred as ITS) from ribosomal RNA (rRNA) gene clusters, among other genes, remain the preferred sequence target for phylogenetic and taxonomic analysis. These two regions, the mitochondrial COI and ITS region, are widely used in all kingdoms, and even they are being used for the generation of a DNA barcode for species identification [34]. Despite this great portfolio of genetic and genomic data, these preferred fragments are only partially available for some phytophagous mites and are almost lacking in their main biological control agent, the phytoseiid mites that seem the ever-lost group (despite recent works in the past years). The use of COI in mites, as well as in other mites of agricultural interest, is conditioned by a limited number of available sequences of this fragment [35], and by the reliability, sometimes questioned, of the sequences published in the databases [16,36,37]. Indeed, until September of 2021, the public repository of sequence data (GenBank) hosts more than one million of sequences from Acari, but less than 10% of these sequences belong to orders containing species of agricultural importance, and less than 1% belong to phytoseiid mites. In this account, only 49 complete mitochondrial genomes are available, being 24 from the superorder Acariformes and 25 from Parasitiformes, and only two of these 25 belong to Phytoseiidae (in this work, we have followed the taxonomic rules as described by Gu et al. [38] and Ruggiero et al. [39], concerning the high taxonomic level of the groups here studied). All these works highlight the great controversy that exists for the identification of Acari species, based on morphological and/ or on molecular data, and that conforms one of the pillars of this study.
In this work, we have used the Acari mitochondrial COI sequence (following the directives of Ros & Breeuwer [37] and Tixier et al. [16]) and ITS sequence to establishing species boundaries and phylogenetic relationships among several clades of Acari, mainly Tetranychidae and Phytoseiidae.

Materials and Methods
2.1. Acari Collection. Table S1 lists the mites (Tetranychidae and Phytoseiidae) collected mainly from Spanish citrus groves, from laboratory rearing colonies or purchased from different commercial biological control suppliers that have been used to obtain ITS (290 specimens from 14 species) and 3' COI (300 specimens from 12 species) sequences, along with those retrieved from GenBank (36 COI sequences from 34 species and 123 ITS sequences from 85 species). Species identities of field-collected individuals were morphologically determined [40,41].
2.2. DNA Extraction, Amplification, and Sequencing. Total DNA was extracted from single individuals following a modified "salting out" protocol as described in previous works [11,42]. DNA pellet was resuspended in 20 μl of LTE and stored at -20°C.
ITS regions and mitochondrial COI fragments were amplified using the primers listed in Table S2. Exact primer pairs used for each species are listed in Tables S3  and S4. Each fragment was amplified in 25 μl reaction containing 1x Taq polymerase buffer (Roche Applied Science, Mannheim, Germany), 200 μM of each dNTP (5 PRIME GmbH, Hamburg, Austria), 2.5 mM of MgCl 2 , 0.5 μM of each primer, 1 unit of DNA Taq polymerase (Roche), and 1.5 μl of DNA template (5-10 ng/μl). Amplifications were performed in Bio-Rad MJ Research Thermal Cycler PTC-100® and consisted of one denaturation step at 94°C for 4 min, 35 cycles at 92°C for 1 min, annealing at 45 or 50°C (depending on the combination of primers) for 1 min, and 72°C for 1 min 30 s, followed by a final extension at 72°C for 10 min. PCR products were run on 2% agarose D-1 low EEO (Pronadisa, Sumilab S.L., Madrid, Spain) gel, stained with ethidium bromide using a molecular weight marker of 50 bp DNA Ladder (Invitrogen, Carlsbad, California, USA), and visualized under UV light. Each PCR product was purified with Illustra™ ExoStar™ (GE Healthcare Life Sciences, Buckinghamshire, UK) prior to sequencing. Between three and ten different individuals of each species were sequenced (SANGER) in both directions using ABI/ PE 3730 DNA Analyzer (Applied Biosystems, Foster City, USA) at the Servei Central de Suport a la Investigació Experimental (SCSIE), Universitat de València (Spain).

Sequence Analyses.
A consensus sequence for each PCR product was obtained using the program STADEN Package [43] and verified belonging to the proper copy of COI or ITS by an initial BlastN or BlastX homology search against NCBI 2 Journal of Zoological Systematics and Evolutionary Research database [44]. As indicated, outgroup sequences were retrieved from GenBank (Table S1). Consensus COI sequences were aligned using GENE-DOC [45], nucleotide sequences were translated using code 5 (invertebrate mitochondrial) and aligned using Blossum62 score table, and set the alignment cost at 20 for constant length, 8 for gap opening, and 4 for gap extension. This amino acid alignment was used to realign mitochondrial COI nucleotide sequences. ITS sequences were aligned using CLUSTALW option implemented in MEGA X [46], with default parameters for DNA weight matrix, gap opening, and extension penalties. Pairwise mean distance matrices for each marker were obtained with MEGA X (Tables S5  and S6). MODELTEST [47], as implemented in MEGA X, was used to determine the best-fit model to each data set, using it as the seed to determine the final gene tree.
Acari COI and ITS alignments were converted to nexus format for Bayesian phylogeny inference with BEAST. Bayesian phylogenies were obtained in BEAST v1.10.4 program [48] with GTR + I + G evolutionary model using a Markov Chain Monte Carlo (MCMC) method, with 1000000 length of chain, 1000 screenlog, and 200 tracelog parameters, Yule speciation process [49], and applying uncorrelated lognormal relaxed clock [50] for COI region and strict clock for ITS marker. Consensus species tree was generated with TreeAnnotator module in BEAST after burning 10% of departing trees with posterior probability limit >0.5. The consensus final species trees were visualized with FigTree v1.4.3, as implemented in BEAST using as outgroup the Aranea species Lycosa coelestis L. Koch (Araneae: Lycosidae) for COI and Pardosa tristis (Thorell), Latrodectus katipo Powell, and Enoplognatha ovata (Clerck) for ITS.

Results
3.1. New Acari Mitochondrial COI Sequences, Alignment, Barcoding, and Phylogeny Inference. We have obtained 29 new Acari mitochondrial COI sequences (Table S3), which were confirmed as Acari COI by >97% homology with blastX against GenBank dataset.
All new sequences, including the outgroup L. coelestis, were unambiguously aligned (taking into consideration the amino acid sequence), not detecting single nucleotide insertion or deletions, nor stop codons indicating that the obtained mitochondrial COI sequences belong to the active mitochondrial copy not to numts. Only three new sequences, all belonging to the same species (Typhlodromus (Anthoseius) rhenanoides Athias-Henriot) (Mesostigmata: Phytoseiidae), were identified as numts and were removed from downstream analysis. On average across all Acari taxa, the AT content was 75%, which is a general feature of arthropods COI [51,52]. Figure 1 shows graphically the variation occurring in the aligned sequence obtained as a fingerprint (Fingerprint software [53]), similar to the barcode of Hebert et al. [34]. From the complete sequence of COI gene, positions 950 to 1450 have been expanded to highlight the heterogeneity in this gene region, corresponding to the focused Acari mitochondrial COI region sequenced in this work.
The sequence alignment containing all Acari COI sequences was used to determine model best fitting to our data, finding that is general time reversible (GTR) with gamma distribution of heterogeneity and invariant ðGTR + G + IÞ. This DNA substitution model was applied in a first round to determine the starting ML tree ( Figure S2), which was lately used to infer the phylogeny with the speciesreduced sequence alignment under Bayesian species tree inference ( Figure 2).
The phylogenetic tree based on Bayesian inference determined that Acari is a monophyletic group, in which Ixodida, Mesostigmata, Prostigmata, and Oribatida were also monophyletic. The first taxon to split from Acari was determined to be Ixodida, followed by Oribatida and groups Prostigmata and Mesostigmata ( Figure 2). Even if this fragment of the COI sequence did not allow ascertaining the depth branching point of Acari, due to low values of posterior probability, this fragment gave enough resolution for species identification based on sequence distance and phylogenetic position. This allowed us to establish a DNA barcode for Acari species identification based on nucleotide distances with known individuals (despite the controversy raised on the use of nucleotide distance for species assignments as discussed below).
3.2. New ITS Sequences, Alignment, and Analysis of Molecular Evolution Patterns. Prior to starting this work, the number of sequences covering the ITS from Acari was below 100 records, belonging mainly to Acariformes (see  The third row (identity) depicts constant, variable, and gaps present in the alignment. The fourth row (heterozygosity) depicts the frequency of nucleotides in each position. And finally, the fifth row (diversity) represents how variable is each site within the sequence by using the Nei and Li [80] nucleotide diversity measure. As can be shown, the region 950-1450 has been expanded, which forms the COI 3' region sequenced in this work. Table S1). Nowadays, with the new 42 ITS sequences obtained in this study (Table S4) and the scientific community effort, the number of available sequences is above one thousand. However, in many cases, these sequences are only partial (covering only ITS1 or ITS2, or partially both covering in full the 5.8S rDNA that separates each ITS). We have used a total of 155 sequences covering the taxonomic groups Ixodida, Mesostigmata, Prostigmata, Astigmata, and Oribatida. The first point that should be highlighted is the extreme difference in fragment size. Almost all Mesostigmata had an average length of 700 bp, except Gamasinae (within Mesostigmata) that had an average of 800-900 bp, the Prostigmata of nearly 1000 bp, the Astigmata nearly 1600 bp, and finally, the Ixodida ranging between 2000 and 6000 bp in length. These extreme length differences would influence the final alignment, as in some groups, ITS1 is shorter than ITS2 and opposite in others. As this fragment is a noncoding region that is not subjected to functional constraints (protein function) other than structural ones, allowing multiple insertion-deletions events. Figure 3 shows the results of fingerprinting as done for the Acari mitochondrial COI region, indicating that ITS alignment is more informative within positions 1100 to 1800 (nucleotide diversity row), due to the large number of Phytoseiidae used in this work. Interspecific uncorrected pairwise distances were determined from the alignment, showing an average interspecific distance of 0:3913 ± 0:0377 within all Acari against the Araneae outgroups used, of 0:5121 ± 0:0959 within Astigmata, of 0:4499 ± 0:2442 within Ixodida, of 0:4389 ± 0:0170 within Oribatida, of 0:1772 ± 0:0844 within Mesostigmata, and 0:4070 ± 0:1440 within Prostigmata (Table  S6). Intraspecific distances were determined in only 9 species, 3 belonging to Prostigmata (P. citri 0:00 ± 0:00 (n = 8); T. evansi 0:0039 ± 0:0019 and 0:0009 ± 0:0016; and T. urticae 0:0096 ± 0:0094 and 0:022 ± 0:044 (for ITS1 and ITS2, respectively)) and another 6 belonging to Mesostigmata (Amblyseius andersoni (Chant) 0:0056 ± 0:0087, E. stipulatus 0:00 ± 0:00 (n = 2), N. barkeri 0:0083 ± 0:0086, Neoseiulus californicus (McGregor) 0:0013 ± 0:0009, Phytoseiulus persimilis Athias-Henriot 0:0038 ± 0:0028, and T. phialatus 0:0072 ± 0:0031) ( Table S6) Figure 2: Evolutionary relationships of Acari taxa based on Cytochrome C oxidase I (COI) sequence alignment. Phylogeny inference was performed using Bayesian analysis (in BEAST) of 41 nucleotide sequences, under GTR + I + G model for DNA substitution, and under Yule speciation process with only one sequence representing each species. Ambiguous positions with less than 50% site coverage were eliminated, rendering only 606 positions in the final dataset. Acari are shown to be a monophyletic group, with all species belonging to Ixodida (in blue), Oribatida (in yellow), Prostigmata (in green), and Mesostigmata (in orange).

Journal of Zoological Systematics and Evolutionary Research
6:93043E-20; and T. phialatus p = 1:4834 E-19 in Mesostigmata), and due to the 100% similarity on several specimens, we were able to reduce the number of sequences to one per species, to obtain the species tree with ITS marker.
As indicated in the materials and methods section, the starting of all sequence alignment (the alignment containing all specimens sequences, including those with 100% similarity) was used to determine the best fitting model, applying the GTR + G + I model in the next steps and establishing  The first row (nucleotides) shows the nucleotide composition for the consensus sequence obtained from the alignment. The second row depicts the heterogeneity of the sequences in each nucleotide position. The third row (identity) depicts constant, variable, and gaps present in the alignment. The fourth row (heterozygosity) depicts the frequency of nucleotides in each position. And finally, the fifth row (diversity) represents how variable is each site within the sequence by using the Nei and Li [80] nucleotide diversity measure. As can be shown, the region from 2600 to 3054 is black (for nucleotides, heterogeneity, and identity) or red (for heterozygosity and diversity) colored, meaning that this region is only present in less than 1% of sequences. Noticeable is the central region of the alignment, which shows the most informative sites. 6 Journal of Zoological Systematics and Evolutionary Research early inference of phylogeny with less robust methods (see Figure S3). The Bayesian phylogeny inference of species tree indicates that Acari form a monophyletic group with a strong support of the branch (Figure 4). Within Acari, only Prostigmata seemed not behaving as monophyletic group, with only two sequences obtained from GenBank belonging to genus Cecidophyopsis Keifer that seems to deserve further study concerning their species status. Oribatida and Astigmata were the early-branching groups in Acari clearly separated by a clade formed by Mesostigmata and with Ixodida as sister group of the main Prostigmata subgroup (Figure 4).

Discussion
The use of DNA sequences has unveiled a great diversity of Acari, even within those described as species synonyms ( [37], and references herein). However, not all people working with Acari are dealing with the same gene or region within the selected gene marker, and the number of whole genomes is still low enough to perform phylogenomic studies covering all the Acari major groups (see [54]). For species identification and delineation, within "Tree of Life" project, it is generally used the "Folmer fragment" of COI gene, which is limited by the primers LCO1490-HCO2198 ( [55]; see Figure 1), or is used also the major subunit of ribosomal RNA (18S rDNA or 16S rDNA, depending on if we consider eukaryotes or prokaryotes). The Folmer fragment is recognized also as a taxonomical character [34], accounting for more than 1.3 million of records in GenBank and has been adopted by some Plant Protection and Quarantine agencies (like USDA-ARS) to control and monitor pest and invasive species [56]. However, the lack of records covering Acari species or not covering their "barcode" region (corresponding to the 3' region of COI gene or ITS) ( [4,16,33,37,57]; and references herein) highlights the need supplied in this work by increasing the number of sequences belonging to agriculture important Acari species, also being able to fit some gaps missing in phylogenies, which is a hot topic in Acarology nowadays [13,54,[58][59][60][61]. However, species delimitation based solely on molecular data is still under controversy, with some groups asking for integrative methods that take into account nucleotide distances, phylogenetic closeness, and morphological evidence, placing efforts to develop analysis tools for this integrative species delimitation [62][63][64][65][66][67][68], none of them being focused on Acari species.
With this study and the Pérez-Sayas et al. [11] work, we have increased the number of sequences covering the Acari mitochondrial COI region and ITS for both Tetranychidae (within Prostigmata) and Phytoseiidae (within Mesostigmata). The sequence analyses of both DNA regions have unveiled their use as a species identification character, with special emphasis on Acari mitochondrial COI in comparison Figure 4: Evolutionary relationships of Acari taxa based on internal transcribed spacer (ITS) sequence alignment. Phylogeny inference was performed using Bayesian analysis (in BEAST) of 96 nucleotide sequences, under GTR + I + G model for DNA substitution, and under Yule speciation process with only one sequence representing each species. Ambiguous positions with less than 50% site coverage were eliminated, rendering only 3003 positions in the final dataset. Acari are shown to be a monophyletic group, as with COI. This phylogenetic inference shows Oribatida (in yellow) as basal group grouping with Astigmata (in purple), whereas Prostigmata (in green) is split into two subgroups, one formed by two sequences that is basal to Mesostigmata (in orange) and Ixodida (in blue), with the main group of Prostigmata being the sister cluster of Ixodida. 200   400  600  800  1057  1   200  400  600  800  1000  1200  1400  1539  1   200  400  600  800  1000  1200  1400  1536  1   1  100  200  300   1   1  50  100  150  200  250  300  350  410   1  50  100  150  200  250  300  350  400  453   1  50  100  150  200  250  300  350  400  453   1  50  100  150  200  250  300  350  400  454   1  50  100  150  200  250  300  350  400  454   1   1  200  400  600  800  1000  1200  1400  1554   50  100  150  200  250  300  350  400  475   1  50  100  150  200  250  300  350  400   1  50  100  150  200  250  300  362   453   1  50  100  150  200  250  300  350  Amblyseius andersoni (Chant) Figure 5: Cytochrome C oxidase I (COI) barcodes of several Acari of agricultural interest. Based on each COI independent sequence, barcodes were obtained and aligned to the approximately nucleotide position with respect to Tetranychus urticae str. red those of Acariformes (with grey background), or aligned with Metaseiulus occidentalis sequence corresponding to Parasitiformes (dark grey lines). Some outgroups have also been included to notice for first view differences. 8 Journal of Zoological Systematics and Evolutionary Research with the Folmer fragment ( Figure 5). This marker showed an interspecific mean distance (0:2760 ± 0:0214) that was orders of magnitude higher than the intraspecific distance, allowing a proper species assignment by comparison with the deposited sequences in GenBank (see Table S5). With the sole exception of E. nicholsi sequences, which either may form a complex of cryptic species within the original area (China and South Asia), or belong to reproductively isolated populations in the limit of speciation [69]. In addition, the lack of differences between populations in some pest species like P. citri would allow to convert these sequences in a proper DNA barcode for species identification purposes (like the ones required Plant Protection agencies controlling borders). However, to properly disentangle species boundaries with the Acari COI fragment sequence, more populations should be analyzed expanding geographical areas and hosts, as mitochondrial inheritance can be distorted by the presence of bacterial symbionts like Wolbachia or Serratia that affect species reproduction [70]. To avoid this drawback of having sex-distorter agents, it is generally accepted the use of nuclear markers like the one used here, the ITS, as an alternative. With the results here provided, increasing the number of sequences of ITS and by determining the intraspecific p-distances, we corroborate its usage as another DNA barcoding region that allows for reliable species identification based on nucleotide p-distances or phylogenetic closeness to known voucher specimens' sequences. These two fragments have shown enough contiguous sequence to develop species-specific primers to set up a PCR marker which is easily implemented in pest management programs, as many of the pest management officers have access to routine low-cost molecular biology laboratories (see [9,11]). In this work, we were able to infer DNA barcodes from ITS sequences, as they show high interspecies and low intraspecific nucleotide distances that altogether with phylogenetic relationships allow the assignment to species (see Table S6). However, species divergence threshold for ITS or COI sequences might be underestimated due to the low number of populations tested, as also occurred in other arthropods groups [71,72]. Our phylogenetic analyses recovered the traditional groupings within the Acari (Figures 2 and 4, and supplementary Figures), including its monophyletic state within the Chelicerata, when using one or several sequences of Aranea or other arthropods (aphids (Hemiptera) as an outgroup of Chelicerata; data not shown). However, we do not recover as monophyletic the superorders Acariformes (Prostigmata + Oribatida + Astigmata) and Parasitiformes (Mesostigmata + Ixodida) as happened in Lozano-Fernandez et al. [59] using phylogenomics. Other authors also recovered this monophyly when using between one to three sequences belonging to the Mesostigmata, or when using Parasitiformes as outgroup of Acariformes ( [35,38,73]; and references herein). In this work, we recover the polyphyly of Acariformes and Parasitiformes with both markers (Figures 2 and 4) similarly to a study that used only nuclear ribosomal genes and secondary structure alignment, despite that these authors only included two Phytoseiidae species within the Parasitiformes [74]. The paraphyletic status of Prostigmata determined with ITS should be taken cautiously; as explained, the phylogenetic signal may be influenced by long-branch attraction, similarly to those obtained for Prostigmata with concatenated mitochondrial proteins [38] or COI Folmer fragment and 18S rDNA in other studies [30,60,75].

Journal of Zoological Systematics and Evolutionary Research
In a depth view, our phylogenetic analysis recovers partially the life traits of Acari, finding in many tree reconstructions, the Oribatida as a basal group, indicating a detritivorous feeding behavior of the first emerging Acari as postulated with morphological-based phylogeny ([2] and references herein) or by molecular dating of Acariformes based on COI Folmer fragment [30].
In conclusion, despite the controversy raised on the use of molecular data to establish species boundaries or phylogenetic position to establish species assignment of unidentified specimens [62-68, 76, 77], the identification of Acari pests and their predatory mites, and the phylogeny of Acari, remains a controversial issue, being forced to rely on sequence data by the reduced number of discriminating morphological characters, and by their similitude with ancient Acari ( [2,78] and references herein; [79]). However, an effort should be made to coordinate the barcoding region used to calibrate both phylogeny and species limits. Both parameters are required for species assignment based on sequence data when the discriminating morphological characters are either difficult to determine or show phenotypic plasticity sometime bigger within populations than between species. As when the success of quarantine measures or pest management protocols relies on the correct identification of pests and natural enemies.

Data Availability
Sequences obtained and used in this study are included in the supplementary information files, and the new sequences obtained are deposited in GenBank (NCBI).

Disclosure
This work is part of C. Pérez-Sayas PhD thesis.

Conflicts of Interest
The author(s) declare(s) that they have no conflicts of interest.

Authors' Contributions
Consuelo Pérez-Sayas and Tatiana Pina contributed equally to this manuscript and should be considered as first authors.

Supplementary Materials
Additional supporting information may be found in the online version of the article at the publisher's website. Table  S1: list of all sequences used in this work. Table S2: list of primers used to generate the 3' COI fragment, ITS1, ITS2, or whole ITS from either Phytoseiidae or Tetranychidae used in this work. Table S3: new Acari mitochondrial COI sequences obtained in this work. Table S4: new ITS sequences obtained in this work. Table S5: estimates of evolutionary divergence distance between Acari COI sequences. This excel file contains several sheets, first one with interspecific distances ± SE, and the others with intraspecific distances ± SE (one sheet per species). Table S6: estimates of evolutionary divergence distance between Acari ITS sequences. This excel file contains several sheets, first one with interspecific distances ± SE, and the others with intraspecific distances ± SE (one sheet per species). Figure S1: evolutionary relationships of Acari taxa based on 43 Acari Cytochrome C oxidase I (COI) fragment sequences inferred by BEAST (Bayesian analysis), including three different sequences from Euseius nicholsi. Figure S2: evolutionary relationships of Acari taxa based on 61 Acari Cytochrome C oxidase I (COI) fragment sequences inferred by MEGA X [46], using the maximum likelihood (ML) method. Figure  S3: evolutionary relationships of Acari taxa based on 157 Acari Internal Transcribed Spacer (ITS) fragment sequences inferred by MEGA X [46], using the maximum likelihood (ML) method. (Supplementary Materials)