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

Authors: Consuelo Pérez-Sayas; Tatiana Pina; Beatriz Sabater-Muñoz; María Antonia Gómez-Martínez; Josep A. Jaques; Mónica A. Hurtado-Ruiz
Journal: Journal of Zoological Systematics and Evolutionary Research
 (2022)

---

## Abstract

Acari harbor numerous minute species of agricultural economic importance, mainly Tetranychidae and Phytoseiidae. Great efforts have been established by means of recovering morphological, molecular, and phylogenetic traits for species identification. Traditional identification still relies on external diagnostic characters, which are limited and usually exhibit large phenotypic plasticity within the species, rendering them useless for species delimitation and identification. We decided to increase the number of sequences of the Acari mitochondrial COI (Cytochrome C oxidase I) marker and ITS nuclear ribosomal DNA region for species identification 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’ identification and phylogenetic studies. Analyses of both DNA regions have unveiled their use as species identification characters, with special emphasis on Acari mitochondrial COI for Tetranychidae and Phytoseiidae species in comparison with theFolmer 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 identification. 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.


---

## Body

## 1. 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 taxon-linked 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–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–18], to the understanding of the dynamic of the pest [19], the establishment of dispersal patterns [19, 20], the determination of population structure [21–27], their food preferences [9, 11], and even to the determination of geographical origin of invasive species producing an outbreak [28–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.


## 2. Materials and Methods

### 2.1. Acari Collection

TableS1 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 TableS2. 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 MgCl2, 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).


### 2.3. 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 database [44]. As indicated, outgroup sequences were retrieved from GenBank (Table S1).Consensus COI sequences were aligned using GENEDOC [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.


## 2.1. Acari Collection

TableS1 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 TableS2. 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 MgCl2, 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).


## 2.3. 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 database [44]. As indicated, outgroup sequences were retrieved from GenBank (Table S1).Consensus COI sequences were aligned using GENEDOC [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.


## 3. Results

### 3.1. New Acari Mitochondrial COI Sequences, Alignment, Barcoding, and Phylogeny Inference

We have obtained 29 new Acari mitochondrial COI sequences (TableS3), which were confirmed as Acari COI by >97% homology with blastX against GenBank dataset.All new sequences, including the outgroupL. 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.Figure 1
Cytochrome C oxidase I (COI) alignment fingerprint. 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 950-1450 has been expanded, which forms the COI 3’ region sequenced in this work.Interspecific uncorrected pairwise distances were determined from the alignments, showing an average interspecific distance of0.2760±0.0214 within all Acari against the Araneae outgroup used, of 0.1823±0.0183 within Ixodida, of 0.2161±0.0585 within Oribatida, of 0.2056±0.0275 within Mesostigmata, and of 0.1492±0.0228 within Prostigmata (Table S5). Intraspecific distances were determined in 8 species, 4 belonging to Prostigmata (with average distance of 0.00±0.00 in Panonychus citri (McGregor) (n=5); Tetranychus evansi Baker and Pritchard 0.0024±0.0034; T. okinawanus Ehara 0.0044±0.0030; and T. urticae Koch 0.0638±0.0193) and another 5 belonging to Mesostigmata (Euseius stipulatus (Athias-Henriot) 0.0132±0.0049; Kampinodromus aberrans (Oudemans) 0.0169±0.0068; Neoseiulus barkeri Hughes and Typhlodromus phialatus Athias-Henriot 0.0139±0.0121) (Table S5). The intraspecific distances were significantly smaller than the interspecific distances of each subgroup (P. citrip=1.95168E‐16; T. evansip=2.20338E‐13; T. urticaep=0.00058 in Prostigmata; E. stipulatusp=1.25981E‐12; and T. phialatusp=4.51045E‐12 in Mesostigmata) allowing us to reduce the number of sequences to one per species, to obtain the species tree with Acari COI marker. The fifth species, E. nicholsi (Ehara and Lee), from Phytoseiidae (Mesostigmata), showed an average p-distance (0.201±0.01) that not was significatively different from the average p-distance in Mesostigmata (E. nicholsip=0.35851), indicating that each of the sequences analyzed may correspond to a cryptic species complex (as shown in Figure S1), or highly isolated populations.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 invariantGTR+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 species-reduced sequence alignment under Bayesian species tree inference (Figure 2).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, underGTR+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).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 (Figure2). 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 TableS1). 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.Figure 3
Internal transcribed spacer (ITS) alignment fingerprint. In this case, the structure of the sequenced fragment has been included at the top of the figure. 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.Interspecific uncorrected pairwise distances were determined from the alignment, showing an average interspecific distance of0.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. citri0.00±0.00 (n=8); T. evansi0.0039±0.0019 and 0.0009±0.0016; and T. urticae0.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. stipulatus0.00±0.00 (n=2), N. barkeri0.0083±0.0086, Neoseiulus californicus (McGregor) 0.0013±0.0009, Phytoseiulus persimilis Athias-Henriot 0.0038±0.0028, and T. phialatus0.0072±0.0031) (Table S6). The intraspecific distances were significantly smaller than the interspecific distances of each subgroup (T. evansip=1.7326E‐14, 1.4578E‐14T. urticaep=1.8441E‐14, 1.14319E‐14 (ITS1 and ITS2, respectively) in Prostigmata; A. andersonip=5.55465E‐20; N. barkerip=6.6676E‐19; P. persimilisp=6.93043E‐20; and T. phialatusp=1.4834E-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 theGTR+G+I model in the next steps and establishing 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).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, underGTR+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.


## 3.1. New Acari Mitochondrial COI Sequences, Alignment, Barcoding, and Phylogeny Inference

We have obtained 29 new Acari mitochondrial COI sequences (TableS3), which were confirmed as Acari COI by >97% homology with blastX against GenBank dataset.All new sequences, including the outgroupL. 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.Figure 1
Cytochrome C oxidase I (COI) alignment fingerprint. 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 950-1450 has been expanded, which forms the COI 3’ region sequenced in this work.Interspecific uncorrected pairwise distances were determined from the alignments, showing an average interspecific distance of0.2760±0.0214 within all Acari against the Araneae outgroup used, of 0.1823±0.0183 within Ixodida, of 0.2161±0.0585 within Oribatida, of 0.2056±0.0275 within Mesostigmata, and of 0.1492±0.0228 within Prostigmata (Table S5). Intraspecific distances were determined in 8 species, 4 belonging to Prostigmata (with average distance of 0.00±0.00 in Panonychus citri (McGregor) (n=5); Tetranychus evansi Baker and Pritchard 0.0024±0.0034; T. okinawanus Ehara 0.0044±0.0030; and T. urticae Koch 0.0638±0.0193) and another 5 belonging to Mesostigmata (Euseius stipulatus (Athias-Henriot) 0.0132±0.0049; Kampinodromus aberrans (Oudemans) 0.0169±0.0068; Neoseiulus barkeri Hughes and Typhlodromus phialatus Athias-Henriot 0.0139±0.0121) (Table S5). The intraspecific distances were significantly smaller than the interspecific distances of each subgroup (P. citrip=1.95168E‐16; T. evansip=2.20338E‐13; T. urticaep=0.00058 in Prostigmata; E. stipulatusp=1.25981E‐12; and T. phialatusp=4.51045E‐12 in Mesostigmata) allowing us to reduce the number of sequences to one per species, to obtain the species tree with Acari COI marker. The fifth species, E. nicholsi (Ehara and Lee), from Phytoseiidae (Mesostigmata), showed an average p-distance (0.201±0.01) that not was significatively different from the average p-distance in Mesostigmata (E. nicholsip=0.35851), indicating that each of the sequences analyzed may correspond to a cryptic species complex (as shown in Figure S1), or highly isolated populations.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 invariantGTR+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 species-reduced sequence alignment under Bayesian species tree inference (Figure 2).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, underGTR+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).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 (Figure2). 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 TableS1). 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.Figure 3
Internal transcribed spacer (ITS) alignment fingerprint. In this case, the structure of the sequenced fragment has been included at the top of the figure. 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.Interspecific uncorrected pairwise distances were determined from the alignment, showing an average interspecific distance of0.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. citri0.00±0.00 (n=8); T. evansi0.0039±0.0019 and 0.0009±0.0016; and T. urticae0.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. stipulatus0.00±0.00 (n=2), N. barkeri0.0083±0.0086, Neoseiulus californicus (McGregor) 0.0013±0.0009, Phytoseiulus persimilis Athias-Henriot 0.0038±0.0028, and T. phialatus0.0072±0.0031) (Table S6). The intraspecific distances were significantly smaller than the interspecific distances of each subgroup (T. evansip=1.7326E‐14, 1.4578E‐14T. urticaep=1.8441E‐14, 1.14319E‐14 (ITS1 and ITS2, respectively) in Prostigmata; A. andersonip=5.55465E‐20; N. barkerip=6.6676E‐19; P. persimilisp=6.93043E‐20; and T. phialatusp=1.4834E-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 theGTR+G+I model in the next steps and establishing 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).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, underGTR+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.


## 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–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–68], none of them being focused on Acari species.With this study and the Pérez-Sayaset 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 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].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 toTetranychus 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.Our phylogenetic analyses recovered the traditional groupings within the Acari (Figures2 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].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.



---

*Source: 5317995-2022-06-21.xml*