Australian Bryobia mites (Trombidiformes: Tetranychidae) form a complex of cryptic taxa with unique climatic niches and insecticide responses.

BACKGROUND: Bryobia (Koch) mites belong to the economically important spider mite family, the Tetranychidae, with >130 species described worldwide. Due to taxonomic difficulties and most species being asexual, species identification relies heavily on genetic markers. Multiple putative Bryobia mite species have been identified attacking pastures and grain crops in Australia. In this study, we collected 79 field populations of Bryobia mites and combined these with 134 populations that were collected previously. We characterised taxonomic variation of mites using 28S rDNA amplicon-based DNA metabarcoding using next-generation sequencing approaches and direct Sanger sequencing. We then undertook species distribution modelling of the main genetic lineages and examined the chemical responses of multiple field populations.
RESULTS: We identified 47 unique haplotypes across all mites sampled that grouped into four distinct genetic lineages. These lineages have different distributions, with three of the four putative lineages showing different climatic envelopes, as inferred from species distribution modelling. Bryobia mite populations also showed different responses to a widely used insecticide (the organophosphate, omethoate), but not to another chemical (the pyrethroid, bifenthrin) when examined using laboratory bioassays.
CONCLUSION: Our findings indicate that cryptic diversity is likely to complicate the formulation of management strategies for Bryobia mites. Although focussed on Australia, this study demonstrates the challenges of studying Bryobia and highlights the importance of further research into this complex group of mites across the world.

INTRODUCTION

Bryobia (Koch) mites belong to the economically important spider mite family, the Tetranychidae. They are widely distributed across North America, South America, Europe, Africa, New Zealand, Asia and Australia, , , , , , with >130 species described worldwide, although these are likely to contain synonyms and species with overlapping descriptions. Species have previously been described based on morphological characters and plant host associations. , , , , However, suitable morphological characteristics are limited in these small mites, making identification difficult, , and species descriptions based primarily on plant‐host associations are often problematic.

To further compound difficulties with taxonomic resolution of Bryobia species, the majority of described species reproduce asexually by thelytokous parthenogenesis, with populations made up entirely of females; males have only been described from a few species. , , Bryobia and other spider mites are ancestrally haplodiploid, where fertilised eggs give rise to diploid females and unfertilised eggs give rise to haploid males. , , Many Bryobia species are also known to be infected with Wolbachia and Cardinium, maternally inherited intracellular bacteria known to alter the reproductive mode of infected hosts to enhance their spread. , , , Wolbachia infections in at least two Bryobia species induce parthenogenesis through a process known as gamete duplication, and it is likely that Wolbachia is responsible for parthenogenetic reproduction within this genus more generally.

Bryobia mites are highly polyphagous, with numerous species being considered horticultural pests attacking a variety of fruits trees, while others feed on various grasses, herbaceous plants and grain crops. , , Mites feed on plant material by puncturing plant cells with their mouth parts and sucking out the contents, leaving etiolated trails composed of whitish‐grey spots. , , In Australia, Bryobia mites are significant pests of winter grain crops and pastures, and their importance as pests in these agricultural commodities has increased over recent decades. , , They are known to damage canola, wheat, oats, lucerne, clover and pulse crops, and tend to be most damaging at the establishment phase of plant development. , This is similar to other parts of the world, including South Africa, India, Europe , , and North America. In Australia, Bryobia consist of a species complex comprised of at least 10 putatively discrete taxa based on morphology and DNA sequence information. , , Seven of these putative species (B. sp. I, B. sp. IV, B. sp. VII, B. sp. VIII, B. sp. IX, B. sp. X and Bryobia praetiosa Koch) are known to be present in Australian pastures and grain crops. Of these, B. sp. VIII, B. sp. IX and B. sp. I are believed to be the most common, although this is based on limited sampling, mostly from the state of Victoria. Little is known about the broader biogeographic distribution and abundances of different Bryobia taxa across Australia.

The presence of cryptic species can complicate pest management decisions, as separate control strategies may be needed for individual species. In Australia, Bryobia mites are largely controlled using insecticides through foliar sprays or insecticide seed treatments, with the organophosphate, omethoate, and the pyrethroid, bifenthrin, being the most commonly used. However, there have been increasing concerns around chemical control difficulties in the field involving Bryobia mites. , Although Bryobia mites can possess inherent tolerances to chemicals, it is unclear how these might differ between species, nor is it known whether different species have experienced varying levels of selection for increased insecticide tolerance. Such issues have previously arisen in cryptic mite species from the genus Penthaleus, which are also pests in Australian grain crops and pastures. , , To date, only the chemical sensitivity of Bryobia sp. I has been investigated in Australia.

In mites, DNA‐based markers widely used to identify intraspecific and interspecific variation include mitochondrial DNA (mtDNA) and nuclear ribosomal DNA (rDNA); the mitochondrial cytochrome c oxidase subunit I gene (COI), and the internal transcribed spacers (ITS1 and ITS 2) are the most commonly used markers. , Sequences of the COI gene and ITS regions have been widely utilized to delineate and identify morphologically cryptic species, particularly within the Tetranychidae. , , , , , , , , These genes along with the D1 region within the nuclear 28S rDNA gene have previously been used to determine the phylogenetic relationship among Bryobia mite species in Europe, Australia and Africa.

In this paper, we investigated the species status and distribution of Bryobia mites in Australian pastures and winter grain crops. We undertook amplicon‐based DNA barcoding with 28S rDNA using next‐generation sequencing (NGS) approaches and Sanger sequencing on Bryobia mites collected from winter grain crops, pastures and roadside vegetation throughout Western Australia (WA), South Australia (SA), Tasmania (Tas), New South Wales (NSW), Queensland (Qld) and Victoria (Vic). Sequence variation was analysed using Bayesian inference and hierarchical clustering of genetic distances, which led us to infer that Australian Bryobia mites form at least four distinct genetic lineages. Using presence‐absence data and spatial climatic data, we then modelled the Australian distributions of these four lineages and identified climatic variables that were correlated with their distributions. We also explored the niche overlap of lineages commonly observed in this study and tested the responses of multiple mite populations to commonly used insecticides to further determine if there are phenotypic differences between Bryobia lineages.

MATERIALS AND METHODS

Mite collections

In this study, we utilized historical collections of Bryobia mites and undertook targeted field surveys in parts of Australia that had not been well sampled previously. One hundred and thirty‐four populations of Bryobia mites collected from Vic, NSW, SA and the southern part of WA were utilised from Arthur et al. Additionally, field collections of Bryobia mites were undertaken between April 2017 and October 2019, mostly targeting previously unsampled regions, including Tas and southern Qld (Fig. 1). Mites were collected along roadsides or within pasture fields and winter grain crops. Sites were targeted based on known plant hosts of Bryobia mites, , , including cereals, lucerne (Medicago sativa), capeweed (Arctotheca calendula), clover (Trifolium spp.), barley grass (Hordeum leporinum), wild radish (Raphanus raphanistrum), wild oats (Avena spp.), vetch (Vicia spp.) and canola (Brassica napus). Mites were collected via suction using a Stihl SH55 blower vacuum with a fine gauze mesh placed over the end of the vacuum tube. Samples were placed into plastic containers and transferred to the laboratory. Bryobia mites were then distinguished from other mite species based on morphology (colour, size and arrangement of dorsal setae) using a stereo microscope at ×20 magnification and stored in 100% ethanol at −20 °C for DNA barcoding. Global Positioning System coordinates and plant host were recorded at each sampling location. Collections undertaken between 2017 and 2019 yielded an additional 79 field populations of Bryobia mites (Fig. 1).

Figure 1

Map of Australia showing the distribution of sample locations for Bryobia mites used in this study. Black circles indicate samples collected in Arthur et al. Grey circles indicate locations where Bryobia mites were collected in this study.

DNA extraction, sequencing and bioinformatics

We used an amplicon NGS metabarcoding approach following Fadrosh et al. to determine Bryobia haplotypes present within the 134 samples collected by Arthur et al. and the 79 populations collected in this study. Approximately 20 mites per location were pooled into a single 2 mL tube, homogenised by hand and DNA extracted using DNeasy Blood & Tissue Kits (Animal Tissue Spin‐Column Protocol; Qiagen, Victoria, Australia). Using the primers described in Ros et al., we amplified a region of the 28S rDNA nuclear region from each DNA extraction sample using a two‐step PCR protocol. The first‐round PCR primers contained a (5′–3′) universal adaptor sequence and the forward and reverse primers were from Ros et al. : nextera_28S_D1F‐TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGACCCSCTGAAYTTAAGCAT and nextera_28S_D1R‐GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGAACTCTCTCMTTCARAGTTC. Reaction matrices contained 2 μL DNA, 1× PCR buffer, 3 mM MgCl2, 0.3 μM each primer and 0.2 U of KAPA Plant DNA Polymerase (Kapa Biosystems, Massachusetts, USA) in a final volume of 10 μL. Reaction conditions were as follows: one cycle at 95 °C for 3 min, 40 cycles at 95 °C for 20 s, 48 °C for 15 s, 72 °C for 30 s, and one cycle at 72 °C for 1 min. Two technical replicates were performed for each extraction sample. Following PCR, replicate reactions were pooled together in equal volumes. PCR products were then cleaned with ExoSAP‐it (Applied Biosystems, Massachusetts, USA).

Second‐round PCR primers contained the Illumina p5 or p7 binding region, an 8 bp index sequence and a universal adaptor sequence. Unique forward and reverse index combinations were used for each population sample. Reactions contained 2 μL of cleaned PCR product, 0.5 μM of each primer and 1× Phusion Hot Start II High‐Fidelity PCR Master Mix (Thermo Fisher Scientific, Massachusetts, USA) in a final volume of 20 μL. Reaction conditions were as follows: one cycle at 98 °C for 1 min, 15 cycles at 98 °C for 10 s and 72 °C for 45 s, and one cycle at 72 °C for 10 min. Samples were pooled and then subjected to size selection to remove nonspecific fragments outside of the desired size range using Sera‐Mag SpeedBead Carboxylate‐Modified Magnetic Particles (GE Healthcare Life Sciences, Massachusetts, USA). Sequencing was performed on an Illumina MiSeq platform (Illumina, California, USA) using 300 bp PE chemistry at the Monash Health Translation Precinct (Victoria, Australia). Twelve controls (PCR negative containing water and no DNA) were included in the library.

For some population samples, only a small number of Bryobia mites (fewer than five) were available. For these samples, DNA extraction (DNeasy Blood & Tissue Kits, as above), PCR amplification and Sanger sequencing (ABI 3730xl, Macrogen, South Korea) were performed on individuals using the 28S_D1_F and 28S_D1_R primers as per Ros et al.

Bioinformatic analyses were performed with a custom analysis pipeline that incorporated the software programs VSEARCH v2.9.0 and cutadapt. Forward and reverse reads were first merged before trimming to remove primer sequences. Reads were then de‐replicated and those with an abundance of less than 10 were excluded from further analyses. Reference sequences for phylogenetic reconstructions were obtained from data deposited by Ros et al. and Arthur et al., which included 28S sequences for Bryobia, with the Petrobia genus as an outgroup (GenBank accessions: GU979863‐72, EU487049‐65).

Sequence analyses

Multiple sequence alignments and construction of a phylogenetic tree were performed using Seaview v4.7, with neighbour‐joining methods to obtain a distance‐based tree and pairwise genetic distances between haplotypes. Preliminary analyses suggested that a small number of samples included traces of DNA from the Petrobia genus, grouping with Petrobia sequences used as outgroups. We therefore proceeded with only those haplotypes that fell within the Bryobia subtree. On this dataset, a further filtering step was conducted to retain only those haplotypes at a frequency ≥5% (representing ~1 mite in a sample of 20). Following this, sequences from the Sanger and NGS datasets were combined and filtered to retain only unique haplotypes. Phylogenetic reconstruction was undertaken using Bayesian inference (BI) in Mr Bayes v3.2. We first determined the optimum nucleotide substitution model in Modeltest 3.7 using the Akaike's Information Criterion. The same BI settings were then implemented as in Arthur et al., with five million generations (convergence indicated by average standard deviation of split frequency values <0.01). After discarding the burn‐in, remaining trees were summarized in a consensus tree with posterior probabilities as nodal support.

Given the large number of unique haplotypes and the relatively low nodal support from the BI analysis, we also undertook a clustering analysis in R. The 28S rDNA gene alignment was imported into R using the ape package. Only Bryobia species identified from previous analyses and unique operational taxonomic units (OTUs) identified in this study were analysed, excluding the outgroup Petrobia. A genetic distance matrix between pairwise sequences was constructed using the dist.dna function from ape, with a K80 mutational model and no‐missing data (no pairwise deletion). The hclust function was used to visualise pairwise genetic distances.

The genetic distance matrix was then subjected to a principal coordinate analysis (PCoA), which decomposes the variance in the distance matrix to derive new orthogonal axes that summarise pairwise relationships among sequences. PCoA was performed using the PCoA function within the ape package. To identify those principal coordinate (PCo) axes that explained significant variation in the genetic distance matrix, relative to random expectations, we implemented a broken stick method, using the bsdimension function from the pcdimension package. We tested the total number of positive eigenvalues. Only the first two eigenvalues were identified as being significant. Therefore, the first two PCo axes (PCo1 and PCo2) were used as variables to identify putative taxonomic clusters. Four distinct clusters were observed, which we herein refer to as lineages 1, 2, 3 and 4. Visualisation of our clustering analyses was performed using the ggplot2 and ggtree packages.

Distribution models

We developed species distributions models utilising presence data for the Bryobia mite lineages identified in this study, in addition to absence data for 282 locations taken from Arthur et al. Species distribution models were fit using the boosted‐regression‐tree (BRT) approach, which is a machine learning technique used for classification problems that produces a prediction from an ensemble (boosting) of many weaker models (decision trees). , We chose a bagging size of 0.5, which is appropriate for smaller data sets and implemented the model in R using the code provided by Elith et al. A k‐fold cross‐validation procedure was used to train (90%) and test (10%) each model, whereby the data was randomly divided into 10 subsets, with each subset sequentially omitted from the training process and used to test the model predictions. A learning rate of 0.001 was used, resulting in models with 1500–3050 trees depending on the Bryobia lineage. The relative contribution of each predictor variable to the model was estimated using standard practice. ,

Predictor variables for our BRT models consisted of a subset of four standard BIOCLIM variables obtained from the WorldClim database, which are derived from average monthly (1950–2000) temperature and precipitation data. The variables describe means, trends and seasonal variations of temperature and precipitation, which are more likely to represent physiological limits for species. To enable comparison of estimated environmental responses, we used the same BIOCLIM variable subset across all models fitted to each Bryobia lineage. This subset was selected by examining pair‐wise Pearson's correlation coefficients of predictor variables across the sampling locations in the study. BIOCLIM variables with the highest correlation were removed stepwise until no pairwise correlation coefficient exceeded 0.5 (preferencing lower BIOCLIM variables). This resulted in four BIOCLIM variables that were used across all lineages: BIO01 – Annual Mean Temperature; BIO3 – Isothermality (Mean Diurnal Range/Temperature Annual Range) (×100); BIO9 – Mean Temperature of Driest Quarter; and BIO12 – Annual Precipitation.

This study focused on the winter cropping regions of Australia, largely excluding the large interior regions of Australia's rangelands. Thus, to avoid model extrapolation, the analysis was restricted to areas outside of Australia's rangelands. This was achieved using Australian rangeland boundaries as defined by the Australian Government's Department of Agriculture, Water and the Environment, and the Australian Collaborative Rangelands Information System (ACRIS), reducing the total number of locations used in our climate suitability analysis from 495 to 454.

Niche overlap

Following Schoener, we calculated the niche overlap between mite lineages and as the proportional overlap of presence observations across all sample sites where at least one lineage was recorded or , where and are the relative frequencies of lineage and respectively, across site .

Chemical responses of mites

To examine the chemical responses of Bryobia mites, five populations were collected from fields and roadsides from Vic and NSW in July 2019 (Table 1). Mites were collected via suction as described above using a Stihl SH55 blower vacuum. Populations were placed into plastic containers with paper towelling and vegetation, and then transported back to the laboratory. Mites were kept cool prior to undertaking insecticide bioassays.

Table 1

Collection details of Bryobia populations collected and screened against omethoate and bifenthrin using laboratory bioassays

Population	Latitude	Longitude	Date collected	Plant host
Manangatang	−35.056	143.003	03 July 2019	Medicago sativa
Coolamon	−34.813	147.145	04 July 2019	Broadleaf grasses (Gramineae)
Cootamundra	−34.443	147.966	04 July 2019	Urtica dioica
Harden	−34.659	148.184	05 July 2019	Urtica dioica
Beggan	−34.636	148.294	05 July 2019	Urtica dioica

Two commonly used insecticides, the organophosphate omethoate (Le‐Mat, Cheminova, NSW, Australia) and the pyrethroid bifenthrin (Talstar 250EC, FMC Australasia, NSW, Australia), were used to generate sensitivity data. Le‐Mat contains 290 g L−1 of omethoate and the solution representing the recommended field rate for Bryobia mites was 348 mg a.i. L−1. Talstar 250 EC contains 250 g L−1 of bifenthrin and the solution representing the recommended field rate for Bryobia mites was 200 mg a.i. L−1 of bifenthrin. Bioassays were performed following the methodology described in Arthur et al. but adapted for plastic vials. Seven to nine concentrations of omethoate ranging from 1 × 10−5 to the field rate and eight concentrations of bifenthrin ranging from 1 × 10−5 to the field rate, with the addition of the surfactant Tween at 0.01% (v/v) (Ecoteric 20; Thermo Fisher Scientific, Victoria, Australia), were serially diluted and tested against each population. For each insecticide concentration, approximately 10 mL of solution was poured into a 15‐mL plastic vial and swirled to ensure complete coating, with excess liquid discarded. Between six and eight replicate vials were coated for each concentration and left to dry upside down overnight. The control vials were treated in the same manner, except water (with the addition of 0.01% Tween) was used. Eight Bryobia adults were then placed into each replicate vial along with two leaves of common vetch (Vicia sativa), which provided food and humidity. Vials were placed at 18 °C for 8 h. After this time, individuals were scored as alive (moving freely), incapacitated (inhibited movement) or dead (no movement over a 5‐s period). After scoring, mites from each population were separately placed into Eppendorf tubes containing 100% ethanol and frozen at −20 °C. Fourteen individual mites from each population were later Sanger sequenced to identify Bryobia haplotypes following the methods described above.

As incapacitated individuals invariably die and therefore do not contribute to the next generation, we pooled incapacitated individuals with dead individuals for data analysis across all bioassays. Mite mortality in each bioassay was modelled as a binomial response variable , where is the probability of mortality and is the number of mites in each replicate. The probability of mortality was related to linear covariates with a logistic link function where is chemical dose (mg a.i. L−1) and and are intercepts and slope coefficients, respectively, for each population . , Variance was robustly estimated by allowing for overdispersion.

Each chemical product was analysed separately, with differences in population effects (intercepts) assessed using an analysis of variance (ANOVA) and then examined in pairwise comparisons using Tukey's honest significant difference (HSD) method at the 0.05 significance level, which corrects for type I errors when performing multiple hypotheses tests. LC50 values (with 95% confidence intervals) and slope coefficient estimates were computed. All analyses were performed using the R language for statistical computing.

RESULTS

Bryobia haplotypes

In total, mites from almost 200 locations across Australia, representing approximately 1388 individual Bryobia, were sequenced for the partial (~355 bp) 28S rDNA region using either the pooled NGS amplicon approach (1260 mites, 167 locations) and/or the individual mite Sanger sequencing approach (128 mites, 29 locations). Population samples sequenced included mites from 62 locations in NSW, 42 locations in WA, 41 locations in Vic, 23 locations in Qld, 14 locations in SA and seven locations in Tas. Phylogenetic analyses included 272 bp of the 28S rDNA region, allowing the inclusion of sequences from Ros et al. and Arthur et al. We obtained 47 unique haplotypes from these sequences after filtering. The phylogeny generated from these sequences and the BI node support are shown in Fig. 2. Of the 47 haplotypes, seven were identical to the putative species identified in Arthur et al. (B. sp. I, B. sp. IV, B. sp. VII, B. sp. VIII, B. sp. IX, B. sp. X and B. praetiosa). Similar support to the topology found in Arthur et al. was also found here, with the same two major clades identified: clade 1 contained putative species B. sp. I, B. sp. IV, B. sp. V, B. sp. VII, B. sp. VIII, B. sp. IX, and B. sp. X and clade 2 contained B. praetiosa, B. kissophila, B. berlesei, B. sarothamni, B. rubioculus, B. sp. III and B. sp. VI. Numerous new haplotypes (indicated here as different OTUs) were also found within or in‐between each of the seven putative species from Arthur et al. With the new OTUs identified in this study, there was less support for some of the putative species identified previously, as indicated by the BI posterior probabilities for nodes (Fig. 2).

Figure 2

Bayesian inference tree for 28S rDNA sequence showing relationships between unique haplotypes (OTUs) from this study and putative species identified in Ros et al. and Arthur et al. as well as other related taxa. The broad genetic lineages identified through clustering analysis in this study are shown. Numbers at nodes indicate posterior probabilities (nodal support). Genbank accession numbers are provided in parenthesis.

The PCoA yielded two significant axes that explained the majority of the variation in pairwise genetic distances among haplotype sequences; PCo1 explained 59.6% of the variation and PCo2 explained 26.7% of the variation. Other PCo axes explained less than 5% of the variation in genetic distances and were not significantly different from random expectations. There were four distinct clusters identifiable in PCo space (Fig. 3), which is largely congruent with the BI tree, but some clusters include separate taxa with strong nodal support (Fig. 2). Cluster 1 (herein referred to as lineage 1) included B. sp. I, B. sp. V, B. sp. IX, B. sp. X and 28 OTUs, while cluster 2 (herein referred to as lineage 2) included B. sp. VII and four OTUs. Cluster 3 (herein referred to as lineage 3) is taxonomically diverse and included B. praetiosa and OTU43. The clade where this cluster is located also included B. kissophila, B. sarothamni, B. cf. rubrioculus, B. sp. III and B. sp. VI, but importantly showed deep nodes with bootstrap values that support the separation of these taxa (Fig. 2). Cluster 4 (herein referred to as lineage 4) included B. sp. IV, B. sp. VIII and seven OTUs.

Figure 3

A hierarchical clustering tree of genetic distances between Bryobia haplotypes (OTUs) identified from Arthur et al., Ros et al. and this study. Inset: scatter plot of sequences in PCo space. The labels of the tree and the points in the scatter plot are coloured with respect to K = 4 groups, as inferred from K‐means analysis of PCo1 and PCo2. Genbank accession numbers are provided in parenthesis.

Our analysis raises doubt over the differentiation of the putative Bryobia species identified in Arthur et al. given some previously defined species clustered together in this study which incorporated many more sequences. Instead, Bryobia mites in Australia exhibit a hierarchically structured species complex of four major clades, with several morphologically identified species and putative species nested within one of these clades (lineage 3). This is further supported by the genetic distances within and between the previously defined species of Bryobia sequenced for the 28S rDNA gene (Table S1). The pairwise distances among B. sp. I, B. sp. V, B. sp. IX and B. sp. X (lineage 1) ranged from 0.7% to 2.6%, while the genetic divergence between other putative species were all >7.1%. Pairwise distances between B. sp. VII (lineage 2) and other putative species were substantial, ranging from 7.1% to 9.1%, while pairwise distances between B. praetiosa (lineage 3) and other putative species ranged from 7.9% to 11.2%. The genetic distances among B. sp. IV and B. sp. VII (lineage 4) were <0.6%, with divergence between other putative species in the range 9.9–16.1% (Table S1).

Distribution of mite lineages and environmental suitability

Distributions of the four Bryobia lineages from sampling locations across Australia are shown in Fig. 4. The most common genetic lineages detected across the 454 locations were lineage 1 (127 locations) and lineage 4 (88 locations). Lineage 2 was found at 15 locations. Mites comprising lineage 3 were only found at five locations. Lineage 1 was found across WA, SA, Vic, NSW, Tas and Qld, while the distribution of lineage 4 was slightly more restricted and found in WA, SA, Vic and NSW. Bryobia lineage 2 was restricted to south‐eastern Australia (Vic, NSW and Tas), while the Bryobia clade comprising lineage 3 was only found in Vic and the southern region of NSW.

Figure 4

Distributions and boosted regression tree (BRT) model outputs of climate suitability for Bryobia mite lineages in Australia at a resolution of 0.05° (~5 km). (a) Bryobia lineage 1, (b) Bryobia lineage 2, (c) Bryobia lineage 3 and (d) 4 Bryobia lineage 4. The shading represents the predicted probability of presence based on BIOCLIM variables and fitted BRT models. Black circles represent locations where each lineage was detected.

Environmental suitability maps based on the species distribution models are also depicted in Fig. 4. With the exception of the Bryobia taxa falling into lineage 3, the species distribution models fit the data well, achieving a receiver operating characteristic (ROC) score of 0.88–0.96 for the training data and 0.73–0.82 for the cross‐validation (cv) dataset depending on the lineage. A suitable distribution model could not be fitted to lineage 3, which is not surprising given the relatively few records detected here for B. praetiosa and related OTUs. BIO12 (annual precipitation) was generally the most influential variable in the models, contributing 11.5–46.9% to the models, while BIO03 (isothermality) contributed the least at 10.9–25.3% (Fig. S1). The low overlap in ecological niche observable in Fig. 4 is also demonstrated by the low values calculated for Schoener's index (Table S2), with the exception of lineages 1 and 4, which exhibited a moderate overlap of 0.35.

The estimated environmental responses varied for mite lineages (Fig. 5). The environmental response of Bryobia lineage 1 tended to increase with higher annual mean temperatures (BIO01) but decrease with higher mean temperature of the driest quarter (BIO09) (Fig. 5). In contrast, the response of lineage 2 tended to increase with annual precipitation (BIO12) and decrease with high annual temperatures (BIO01), leading to a predicted distribution concentrated in cooler areas on the Australian mainland and most of Tas (Figs 4 and 5). Similar to lineage 1, Bryobia lineage 4 was found across a wide range, with the exception of eastern coastal regions, which appears to be driven by the strong negative relationship with annual precipitation (BIO12) (Figs 4 and 5).

Figure 5

Partial dependency plots for four BIOCLIM predictors for Bryobia lineages 1, 2 and 4, as determined from a boosted regression tree model (other variables fixed at their average values). Lines are smoothed using a generalised additive model regression spline. The y axes are on a logit scale.

Chemical responses of Bryobia mites

Control mortality of mites was less than 2% across all insecticide laboratory bioassays (Fig. 6). LC50 values for bifenthrin ranged from 5.59 [95% confidence interval (CI) 0.37–84.68] to 17.12 (0.01–22 668.9) mg a.i. L−1 and regression coefficients from 0.475 (±0.061) to 2.668 (±3.538). LC50 values for omethoate ranged from 0.12 (0.03–0.61) to 4.87 (1.65–14.31) mg a.i. L−1 and regression coefficients from 1.058 (±0.599) to 1.735 (±0.960) (Table 2).

Figure 6

Dose–response curves for field populations of Bryobia mites after exposure to (a) bifenthrin and (b) omethoate for 8 h in laboratory bioassays. Error bars represent standard error of the mean. Manangatang and Coolamon are putative species B. sp. IX (i.e. lineage 1). Cootamundra, Harden and Beggan are putative species B. sp. VIII (i.e. lineage 4).

Table 2

LC50 values [and 95% confidence intervals (CIs)] and regression coefficients (and standard errors) for Bryobia mite populations when exposed to omethoate and bifenthrin for 8 h

Insecticide	Population	Bryobia lineage	LC50 value (95% CI) mg a.i. L−1	Regression coefficient b ± SE
Bifenthrin	Coolamon	1	17.12 (0.01–22 668.9)	0.475 ± 0.608
	Manangatang	1	8.15 (0.87–76.03)	2.118 ± 3.078
	Beggan	4	14.68 (2.96–72.80)	2.668 ± 3.538
	Cootamundra	4	9.94 (0.34–290.11)	1.048 ± 1.242
	Harden	4	5.59 (0.37–84.68)	1.115 ± 1.131
Omethoate	Coolamon	1	1.09 (0.20–5.79)	1.275 ± 0.840
	Manangatang	1	4.87 (1.65–14.31)	1.578 ± 0.816
	Beggan	4	0.39 (0.13–1.18)	1.370 ± 0.631
	Cootamundra	4	0.12 (0.03–0.61)	1.058 ± 0.599
	Harden	4	0.16 (0.05–0.45)	1.753 ± 0.960

The genetic lineage of each population, as determined by 28S rDNA sequencing, is also shown.

There was no significant difference in the response of the five mite populations to bifenthrin (χ 2 = 1.75, df = 4, P = 0.78) and there were also no differences when comparing LC50 values which had overlapping 95% CIs (Table 2). This indicates similar bifenthrin sensitivity across all tested populations. Unique regression slopes for each population exposed to bifenthrin did not significantly improve the model (χ 2 = 1.2, df = 4, P = 0.88). In contrast, differences in dose responses were detected between populations exposed to omethoate (χ 2 = 10.08, df = 4, P = 0.04), although unique regression slopes did not significantly improve the model (χ 2 = 0.50, df = 4, P = 0.97). Population differences are further supported by nonoverlapping 95% CIs of LC50 values (Table 2). Manangatang exhibited the highest LC50 value for omethoate, which was 40.5‐fold greater than the LC50 value of the most susceptible population from Cootamundra, while Coolamon had an LC50 value approximately 9‐fold greater than Cootamundra. Mites from Beggan, Harden and Cootamundra had similar dose–response curves and LC50 values for omethoate (Fig. 6 and Table 2).

Sequencing of the 28S rDNA region revealed that Manangatang and Coolamon consisted of mites from putative species B. sp. IX (i.e. lineage 1), while Cootamundra, Harden and Beggan consisted of mites from putative species B. sp. VIII (i.e. lineage 4). These results indicate that mite lineages differ in chemical responses.

DISCUSSION

Bryobia mites are a significant pest of pastures and winter grain crops in Australia and elsewhere, , , , , yet studies on this mite complex are limited. Arthur et al. were the first to apply genetic markers to explore the biology of Australian Bryobia mites and suggested the presence of cryptic species. In this study, we explored the 28S rDNA gene region, which has been successfully used previously to infer species relationships among Bryobia species. Not only did we find the same seven putative species (B. sp. I, B. sp. IV, B. sp. VII, B. sp. VIII, B. sp. IX, B. sp. X and B. praetiosa) identified in Arthur et al., but we also identified a further 40 unique haplotypes compared with only 10 identified by Arthur et al. for the 28S rDNA region. Consequently, there is less support for the clear differentiation of these seven putative species. Rather, our more comprehensive dataset supports the existence of several broad genetic lineages of Bryobia mites within Australian grain crops and pastures. One of these, lineage 3, shows strong nodal support that indicates the presence of multiple species nested within this particular clade. Each of the remaining three lineages may very well represent three distinct species, although delineating species based on genetic data alone (even genome wide analysis) is challenging (see Galtier ), especially in taxa that undergo obligate or facultative asexual reproduction, like many Bryobia.

Bryobia mites are widespread in regions of Australia exhibiting Mediterranean‐type climates (WA, SA, Vic, Tas and NSW), which are characterised by hot dry summers and cool moist winters. However, the Bryobia lineages characterised in this study had different climatic niches based on our species distribution models. Bryobia lineage 1 showed the broadest distribution, being found in all Australian states sampled; it was also the most abundant lineage sampled and had the highest haplotypic diversity (32 haplotypes, 127 locations). Bryobia lineage 2, consisting of five haplotypes, was only detected at 15 sites, and these were restricted to the south‐eastern region of Australia. Mites comprising Bryobia lineage 3 were only found at five sites in Vic and southern NSW, and likely represent the cosmopolitan species B. praetiosa. Similar to lineage 1, Bryobia lineage 4 was abundant among our sample sites, with a widespread distribution across WA, SA, Vic and NSW, although it was absent from eastern coastal regions. The species distribution models highlight that Bryobia mite lineages respond to climatic factors in quite different ways. Lineages 1 and 4 generally expand into Australia's interior rangelands, where precipitation is low. Conversely, lineage 2 appears to be more limited by average temperatures and is predominantly found in cooler areas. The unique ecological niches of these lineages were also supported by the generally low values estimated for Schoener's index. These different environmental responses are likely linked to inherent variation in traits like tolerance of thermal extremes and dry conditions. Such traits have been previously associated with distributions and adaptive shifts in other crop mites in Australia and may be useful in generating hypotheses for future work around interspecific and interclonal variation.

Apart from lineage 3, which includes different putative species, the high degree of haplotype diversity detected in our study in some lineages likely reflects clonal diversity given the majority of Bryobia species reproduce asexually, with males in populations being rare. , Lineage 3 contains 10 haplotypes, most of which were not found in this study. However, the plant hosts for some of these putative species (or haplotypes) were not sampled here. To understand the extent of variation in this lineage and delineate species more accurately, further sampling should focus on other hosts, as well as other geographic locations. For instance, congruous with Arthur et al., we did not detect B. rubrioculus or B. kissophila, which have previously been morphologically identified from Australia. , This likely reflects our sampling approach, which focused on grain crops, pastures and roadside vegetation and did not cover host plants known to be important for these species, such as fruit trees (B. rubrioculus) and English ivy (B. kissophila). , , Also similar to Arthur et al., we did not detect B. graminum, which is surprising given this species has been identified from Australia , and is reported to feed on pasture plants. , Interestingly, Bryobia praetiosa (including OTU‐43) was identified in our study, but only at five locations, which is inconsistent with the notion that this species is a major pest of pasture plants and grain crops in southern Australia. ,

The amplicon NGS metabarcoding approach adopted here for the 28S rDNA region was an efficient method for screening a large number of samples collected from over 167 populations. The method identified 47 haplotypes from >1200 mites, highlighting its efficiency compared with regular Sanger sequencing. However, this approach may under‐represent the true genetic diversity found in Australian Bryobia mites, as the pooled amplicon‐based NGS method and conservative filtering approach adopted here to account for sequencing error will miss some low frequency haplotypes. Nevertheless, our results should reflect the dominant lineages and haplotypes of Bryobia found in Australian grain and pasture environments.

The discovery of multiple lineages of Bryobia mites inhabiting vast areas of agricultural landscapes in Australia signifies possible management challenges for farmers. A series of studies undertaken in the late 1990s and early 2000s showed considerable differences between cryptic Penthaleus mites, which share many life‐history attributes with Bryobia species, including asexuality. The three Penthaleus pest species present in Australia differ markedly in their distributions, , host plants and summer diapause strategies. These species also differ in their response to insecticides, , which is correlated with chemical control failures in the field. , The insecticide bioassays conducted in this study showed significant population differences also exist in Bryobia. Of the five Bryobia populations tested in our study, three were found to belong to lineage 4 and two to lineage 1. The lineage 1 populations had significantly higher LC50 values for omethoate (Manangatang 4.87 mg a.i. L−1, Coolamon 1.09 mg a.i. L−1) compared with all three populations belonging to lineage 4 (which ranged from 0.12 to 0.39 mg a.i. L−1), indicating a higher level of tolerance to this chemical. These results suggest important insecticide tolerance differences between genetic lineages of Bryobia mites, potentially explaining inconsistencies experienced by Australian farmers when attempting to control Bryobia mites with insecticides registered in grain crops and pastures. , Differences in the climatic envelopes inhabited by Bryobia lineages may help inform management strategies in the absence of rapid field diagnostics. For example, lineage 1 is far more common in eastern Australia (NSW and Qld) compared with lineage 4; if the chemical tolerances observed in the present study hold across their observed range, it may be prudent for farmers in eastern NSW and Qld to minimise usage of organophosphates for Bryobia mite control.

In conclusion, this study supports the presence of a diverse species complex of Bryobia mites in Australia. However, contrary to an earlier delineation of seven putative species in Australia, 28S rDNA barcoding of Bryobia populations collected across much of the Australian grain‐growing region supports the differentiation of four broad lineages, one of which (lineage 3) contains a number of distinct species. Additionally, our findings demonstrate that important ecological differences exist between these lineages, including their distribution, climatic suitability and response to insecticides. This highlights the importance of developing management strategies that consider genetic lineages separately, which necessitates the development of tools to readily distinguish Bryobia mites in the field. Currently, this requires considerable skill and experience, and is heavily reliant on genetic markers, with morphological identification by itself unlikely to be satisfactory for separating some lineages. , , , Although focused on Australia, this study demonstrates the challenges of studying Bryobia and highlights the importance of further research into this complex group of mites across the world.

Table S1. Genetic distance (Jukes–Cantor) within and between putative species of Bryobia mites sequenced for the 28S rDNA gene. The putative species identified in Arthur et al.31 are grouped according to their broad genetic lineage identified through clustering analysis in this study. Unique haplotypes of Bryobia from Genbank as well as B. sp. III and B. sp. VI are excluded because of the deep nodes of the phylogenetic tree with bootstrap values that support the separation of these taxa (see Fig. 2)

Table S2. Niche overlap values of the four Bryobia lineages estimated according to Schoener63

Figure S1. The relative contributions to model performance (%) of variables tested as predictors for the distribution of Bryobia lineages 1, 2 and 4 using a boosted‐regression‐tree model.

Click here for additional data file.