The newcomer takes it all: the invader Texas citrus mite, Eutetranychus banksi (Acari: Tetranychidae), displaces the resident relatives in citrus agrosystems

Many studies have emphasized the importance of interspecific competition in shaping natural ecosystem communities. In contrast, few investigations have explored the role of competition in agricultural environments after the arrival of an invasive pest species. We evaluated the ecological impact produced by the invasive Texas citrus mite, Eutetranychus banksi (McGregor), on the resident spider mites Panonychus citri (McGregor) and Eutetranychus orientalis (Klein) on the main citrus crop area in Eastern Spain. Since its arrival in 2013, E. banksi has become the most common and abundant spider mite on citrus, apparently leading to a reduction in the presence and geographic range of the other related species. Competitive relationships were detected between E. banksi–E. orientalis and E. banksi–P. citri pairs using co-occurrence analysis. Furthermore, generalized linear model analysis showed that the probability of finding E. orientalis or P. citri decreases with increasing E. banksi density and vice versa. Principal component analysis and permutational multivariate analysis of variance found competition between these two pairs, and also between the E. orientalis–P. citri pair. Redundancy and variation partitioning analysis revealed how the geographic distribution of the three spider mites is not caused by the environmental conditions, but it is strongly influenced by their colonization history and competitive relationships, since the areas with the highest density of the three species are related to their place of first detection, and do not coincide geographically. Finally, the mechanisms that may be involved in the competitive displacement and the possible future scenarios are discussed.


Introduction
Agricultural ecosystems are usually poor in species richness due to the loss of natural biodiversity, showing simple trophic networks (Duru et al. 2015). Furthermore, these environments are also subjected to constant disturbances through agricultural practices, making them more prone to the establishment of exotic species (Horvitz et al. 1998). Moreover, the use of fertilizers, phytohormones and biostimulants increases food production, but at the same time, enhances the level of available resources in the ecosystem, while a low level of available resources may prevent invasions by alien species (Tilman 1997(Tilman , 1999Tilman et al. 2002).
Phytophagous mites have a high potential to become an invasive species, because their small size make them difficult to detect under routine inspections (Navia et al. 2010). Among them, spider mites, belonging to the family Tetranychidae, cause significant damage in many crops worldwide, and are the most economically important phytophagous mites (Jeppson et al. 1975). In citrus crops these mites feed on the cells of the epidermis and parenchyma of leaves and fruits, reducing the photosynthetic capacity of leaves, as well as affecting the coloration of the fruits, thus producing cosmetic damage that greatly reduces their commercial value (Vacante 2010). Four spider mite species currently occur in the Spanish citrus crop area, the two-spotted spider mite Tetranychus urticae Koch, the citrus red mite Panonychus citri (McGregor), the Oriental red mite Eutetranychus orientalis (Klein), and the Texas citrus mite Eutetranychus banksi (McGregor). Only T. urticae is thought to be native to the Mediterranean basin, while the others invaded the region in the past, the most recent being E. banksi which arrived in 2013 to the Valencia region, the most important Spanish citrus crop area (Ferragut 2016). The behaviour and spatial distribution of these species on the leaves is different. While T. urticae lives on the underside of leaves forming dense colonies, the remaining species prefer to live at lower densities on the upperside, resulting much more likely in competition between Panonychus and both Eutetranychus, as they share the same habitat and compete for space, shelter, and food (Ferragut et al. 2013b).
The Asian species P. citri was the first invasive spider mite to infest the Spanish citrus in 1981 (García-Mari and Del Rivero 1981), spreading rapidly through the Valencian region and becoming a key pest (García-Marí et al. 1983). It was not until 2001 when E. banksi and E. orientalis reached the citrus area in southern Spain (García et al. 2003). Although in both cases the pathway is unknown, E. orientalis is a common species in many countries in Asia and North Africa, while E. banksi is native to the Americas, from southern United States to northern Argentina. Since its arrival in Spain, E. orientalis spread rapidly westward and eastward, reaching the Valencian region in a few years, while the relative E. banksi remained, for about a decade, in a small geographical area. However, in 2013 E. banksi was suddenly detected in citrus orchards in the Valencian southern region, where it most likely arrived through the trade of fruit or plant material (Ferragut 2016). The mite spread rapidly in the region, reaching high populations, producing severe damage to leaves and fruits, and forcing farmers to apply acaricides to maintain populations at low levels.
Since the arrival of Texas citrus mite (TCM) to the Valencian region, farmers and technicians have realized that the spread of the new species coincided with a reduction in the frequency and abundance of P. citri. There is evidence that competition between closely related spider mite species, sharing the same habitat, may limit their geographic range (Fujimoto et al. 1996;Takafuji et al. 1997). Based on field observations, Childers (1992) reported that competition took place between TCM and P. citri in Florida citrus, where TCM has become the dominant species in the orchards since 1955, displacing P. citri.
Competition is a powerful force that regulates the species abundance within communities, but despite their economic importance very few studies have been carried out on the competition between mite species (Reitz and Trumble 2002). Our study aims to evaluate the competitive relationships among the three citrus pests, P. citri, E. orientalis, and E. banksi in the Valencian region, assessing the impact of the recent arrival of TCM in spider mite communities. For this purpose, we address the following objectives: (1) study species co-occurrence using three different approaches: contingency tables, null model, and probabilistic model analysis; (2) define the current spatial distribution of the three mite species using kernel density estimation (KDE); (3) determine the competitive relationships among species using generalized linear models (GLMs), principal component analysis (PCA) and permutational multivariate analysis of variance (PERMANOVA); and finally (4) analyse the relative importance of environmental and spatial factors structuring the community using redundancy (RDA) and variation partitioning analysis (VPA).

Sampling sites
A total of 275 unsprayed citrus orchards were sampled once in a year throughout the three provinces 3173 The newcomer takes it all: the invader Texas citrus mite, Eutetranychus banksi (Acari:… (Castellón, Valencia and Alicante) of the Valencian region from September to December, over a 4-year period (2017)(2018)(2019)(2020). This sampling period was chosen since previous studies have shown that the three mite species have coinciding population peaks in these months with little differences between regions due to climatic differences (Ferragut et al. 1988;Ledesma et al. 2011;Monzó et al. 2016). During the first 2 years, 153 orchards with symptomatic leaves and 95 orchards without symptoms were sampled. In the last 2 years TCM moved northwards, entering the province of Castellón; therefore, during that period 27 additional orchards were sampled in that province ( Fig. 1 in Supplementary Material 1). We used a chi-square goodness of fit test (Gardener 2012) to assess whether the number of samples per province (sampling effort) was significantly different from the expected value according to the citrus area per province.
We collected 20 leaves per orchard that exhibited signs of spider mite damage to determine whether each of our three species were present or absent. This sample size was considered sufficient for detecting the presence/absence of spider mites based on prior surveys of tetranychid mites in the same area (Abad-Moyano et al. 2009). Leaves were placed in plastic bags and transported to the laboratory inside a portable cooling box. In the laboratory, the occurrence of P. citri and Eutetranychus spp. was confirmed under a stereomicroscope. Eutetranychus banksi and E. orientalis cannot be differentiated under stereoscopic microscope, and according to Ferragut et al. (2013b), females must be processed, and slide mounted to distinguish them using a compound microscopic. Reliable characters to identify females of both species include the pattern formed by the setae e1 and f1 on the dorsal surface and the number of ventral setae on the coxa of leg II, and is not necessary to inspect the males to differentiate the two species with certainty (Ferragut et al. 2013b). Eutetranychus density on citrus usually varies between 0-10 females/leaf. Thus, 20 females per sample were randomly collected from the leaves and identified to species level in order to obtain presence/absence data for both species.

Occurrence analysis
Three types of analysis were used to assess co-occurrence patterns in the mite community to check for competitive interactions. First, we mapped the occurrence (i.e., absence/presence) of the three mite species (P. citri, E. banksi and E. orientalis) in all 275 orchards, and then we carried out a (1) null model analysis (Gotelli 2000) by computing two binary (presence/absence) indices for the observed community: checkerboards score (C-score) and variance ratio (V-ratio). C-score quantifies the average amount of co-occurrence among all unique pairs of species in the assemblage (Stone and Roberts 1990), and the V-ratio is calculated as the ratio of the variance of the column sums to the sum of the row variances. For a presence-absence matrix, this is the ratio of the variance in species richness to the sum of the variance in species occurrence (Schluter 1984). Once the indices were obtained for the data-matrix, the values were compared with a distribution of index values obtained from 9999 randomly generated null assembled communities using a randomised algorithm with the R-package EcosimR (Gotelli et al. 2015). To generate the null communities, we used two sequential-swap algorithms with the best statistical properties in terms of Type I and II errors: SIM 2 and SIM 10. SIM 2 algorithm randomizes the occurrence of each species among the sites, assuming the sites are equiprobable. It corresponds to a simple model of community assembly in which species colonize sites independently of one another (Gotelli 2000). SIM 10 algorithm randomizes a binary species-data matrix by reshuffling all of its elements equiprobably (Gotelli et al. 2010). Statistical significance of the observed matrix was calculated as the frequency of simulated matrices that had indices which were equal to or more extreme than the observed index, using the frequency α = 0.05 as the significance level in two-tailed tests. In this way, for the C-Score index competitivelystructured communities should have an index significantly larger than expected by chance, located in the upper tail of the randomly simulated communities index distribution. On the other hand, guild communities will have an index smaller than expected by chance and located in the lower tail. However, competitive communities should have a significantly smaller V-ratio index, and guild communities a larger one, than expected by chance. Finally, the standardized effect size (SES) was calculated following the same methodology as in meta-analysis. This metric was calculated as: observed index − mean (simulated indices)/standard deviation (simulated indices).
It scales the results in units of standard deviations, which allows for meaningful comparisons among different studies (Gotelli and McCabe 2002).
However, species do not occur randomly in nature; they all respond to environmental variation, and algorithms for randomizing species presence-absence data have different statistical properties, computational complexities, and biological realism. Therefore, we calculated species frequency and overlap using the R-package eulerr (Larsson 2021), and used two more methods to assess the co-occurrence patterns. In the first one, we used (2) contingency tables and a Chisquare test (Gardener 2012) to assess the competitive interactions between mite species by testing whether the observed frequencies were statistically different from those expected under the null hypothesis of independence. The second one was made using the R-package coocur (Griffith et al. 2016), computing a (3) probabilistic model of species co-occurrence, which uses basic probability theory to measure cooccurrence as the number of sampling sites where two species co-occur in the community data matrix. Through this method we employed combinatorics to determine the probability that the observed frequency of co-occurrence is significantly large and greater than expected (positive association, i.e. mutualism), significantly small and lesser than expected (negative association, i.e. competition), or not significantly different and approximately equal to expected (random association). Within any species matrix, there will be a mixture of competitive, mutualistic, and random species pairs, and in the above null model analysis, all of which contribute to the observed index value. In contrast to this type of analysis, a probabilistic model analyses the co-occurrence between pairs of species instead of the entire community, which is strictly analytical and requires no randomization being a more parsimonious approach (Veech 2012).

Species distribution and GLM competition models
Geographical position at all 275 orchards was registered in datum WGS84 and reprojected to Cartesian coordinates in ETRS89/UTM zone 30N with the R-packages sp (Pebesma and Bivand 2005;Bivand et al. 2013) and maptools (Bivand and Lewin-Koh 2021). Using these sampled locations, kernel density estimation (KDE) based on a Gaussian function was computed to assess the spatial distribution of P. citri, E. orientalis and E. banksi (Baddeley et al. 2015). We estimated the optimal bandwidth using likelihood cross-validation and a Gaussian function to assess the probability density of the three species per square kilometre, in raster grid cells i.e., raster maps with 3284 × 1890 pixel array and 0.1 square km area per pixel. Once the raster maps were obtained, we extracted the density values for each of the sampled locations by interpolation using the spatstat R-package (Baddeley et al. 2015), and using the density values for each mite species, in each sampled orchard, we obtained the data species matrix for the competition models and multivariate analysis. We performed a Mann-Whitney U-test (Mann and Whitney 1947) to compare the mean density value of each mite species with the presence/absence of the remaining species. Finally, we used binomial (logit-link) generalized linear models (GLMs) (R Core Team 2019) to assess the influence of density on the occurrence, analysing the competitive relationships between the three mite species across the 275 citrus orchards sampled.
Environmental and spatial variables Temperature, relative humidity, and rainfall have been shown to be important factors affecting the development, survival, and reproduction of spider mites, while wind has been described as the main mode of long-distance dispersal (van de Vrie et al. 1972). To examine the influence of these environmental variables on the distribution and competitive relationships between species, a total of sixteen parameters related to temperature, relative humidity, wind speed, and precipitation were considered ( Table 1 in Supplementary Material 2). The data were obtained from 262 weather stations from the Valencian Meteorology Association (AVAMET) scattered around the Valencian region. The span of years for which data were recorded varies across stations, and, for this reason a common time series was used (2013-2020). The stations record environmental data daily, and the annual average in each station for each variable was calculated ( Fig. 1 in Supplementary Material 2). The geographical position of the 262 climatic stations in datum WGS84 was reprojected to Cartesian coordinates in ETRS89/UTM zone 30N as on the species distribution maps. The location and value for the climatic variables in each station were used to perform a spatial smoothing by kernel weighting (KSW), using the Gaussian function with least-squares cross-validation selected bandwidth. We computed the KSW to calculate values for output cells in a raster grid, as in the previous point ( Fig. 3 in Supplementary Material 2). Finally, values for each sampled location were extracted from climate raster maps by interpolation using the spatstat R-package (Baddeley et al. 2015). In addition to climatic variables, other factors such as the citrus crop density can provide information about species competition. We used QGIS to calculate the centroid in datum WGS84 for each municipality (polygon) of a polygons layer (QGIS DT 2009) ( Fig. 2 in Supplementary Material 2), and they were later reprojected in the same way as the climatic stations. The dataset of citrus crop area by municipality was obtained from CADRECTE (2019), and the value for each municipality was associated with its centroid. The centroid geographical locations and its citrus crop area values were used to perform kernel density weighted estimation (KDWE), to assess the probability density of this variable per square kilometre, using the Gaussian function with the crop area value to weigh the points, and the Scott's rule of thumb to determine the smoothing bandwidth (Baddeley et al. 2015). The KDWE allowed us to obtain a raster map with the same number of cells and area per pixel as the climate maps ( Fig. 3 in Supplementary Material 2), and then values for each sampled location were extracted from the raster map according to the same approach. In this way, it was possible to obtain a dataset containing the values of 16 climatic variables and the variable related to the crop area for each of the sampled citrus orchards, building the environmental data matrix for the multivariate analysis. In addition, considering only the localities where the species were present three more environmental data matrices corresponding to each of the mite species were obtained.
We calculated a set of spatial variables as Moran eigenvector maps (MEMs) from the cartesian coordinates of the 275 sampled orchards with at least one mite species, as well as three datasets, one for each mite species, considering only the localities where each was present (Dray et al. , 2012. For this purpose, four connection networks were built following distance based spatial weighting matrix (SWM) using the minimum connection distance as a threshold to define the connected site pairs (Borcard et al. 2011), and the best choice for the weighting functions was done according to the adjusted p-value (Bauman et al. 2018). These analyses were performed with the R-package adespatial (Dray et al. 2021).

Multivariate analysis
To determine the relationship between the three mite species in all 275 sampled locations, the species density data matrix was transformed with the Wisconsin double standardization. This procedure was conducted to normalize sites to the percentage abundance of species, first standardizing species by their maximum value (column maximum), and then sites by their total value (row total) (Legendre and Gallagher 2001). After that, each locality was assigned to a mite species according to its presence. When more than one species was present in the same locality it was assigned to the species with the highest density value, building 3 groups of localities: E. banksi, E. orientalis and Panonychus. Sampled locations with no presence of any of the three mite species were not considered. To test if the grouping was statistically significant, we used non-parametric permutational multivariate analysis of variance (PERMANOVA) (Anderson 2001), considering Chi-square dissimilarity as a measure of distance. The statistical significance was assessed via comparison with 9999 randomized datasets. In addition, analysis of multivariate homogeneity of group dispersions (PERMDISP2) (Anderson 2006) was done to test differential multivariate dispersion (variance) between groups that can affect the results of PERMANOVA. Finally, pairwise PERMANOVA with Bonferroni p-value correction was also carried out to analyse the differences between groups two by two (Martinez-Arbizu 2017). Later, species distribution was analysed with principal component analysis (PCA) showing the significant grouping obtained by the PERMANOVA analysis.
The relative role of the environment and space on the density distribution of the mite species was determined with redundancy analysis (RDA) and variation partitioning analysis (VPA) (Borcard et al. 1992;Peres-Neto et al. 2006), corrected by Moran spectral randomisation (Wagner and Dray 2015). RDA allows us to determine the percentage of variance in a community of species across a set of sites that is due to the variance of a group of variables across those sites. To consider two or more groups of variables, e.g., spatial and environmental variables, the VPA allows us to define what percentage of the variance in the community is due to each of the groups (Borcard et al. 1992). This procedure was applied across the 275 sampled locations with at least one mite species for the whole community, and for each mite species only with locations where the species was present. Multicollinearity among environmental predictors was controlled by the variance inflation factor (VIFs) using VIF < 10 as a preselected threshold (Zuur et al. 2010). Furthermore, we used a forward selection (FS) procedure with two stopping criteria to select environmental and spatial variables in each case (Blanchet et al. 2008) (1) the variance explained by each variable selected had to be significant (p < 0.05), and (2) the procedure stops when the adjusted R 2 accumulated by the variables selected exceeds adjusted R 2 of all the explanatory variables considered. We performed PCA, RDA, PERMANOVA, PERMDISP2 and VIF with the vegan R-package (Oksanen et al. 2020).The R-package adespatial (Dray et al. 2021) was used to perform FS, VPA and Moran spectral randomisation correction.

Sampling effort
Considering the data for 2017 and 2018, the province with the highest sampling effort was Valencia with 160 (64.52%) orchards sampled, followed by Castellón with 47 (18.95%), and Alicante with 41 (16.53%). There were no significant differences between the number of samples per province and the expected value according to the citrus area per province in Valencia (58.23%), Castellón (21.46%) and Alicante (20.32%) (χ 2 = 4.15, df = 2, p value = 0.12). During 2019 and 2020 the sampling effort was increased only in the province of Castellón with 27 orchards, with the aim of updating the geographical distribution of the three mite species (Fig. 1 in Supplementary Material 1).
Patterns of co-occurrence in the mite community were analysed using three different approaches. First, for the null model analysis, C-Score index results under the sim2 and sim10 algorithms showed a nonrandom pattern of species association. Using both algorithms, the C-score observed value (4645.3; p value < 0.05) was significantly higher than the simulated index expected by chance, indicating a negative covariance between mite species, and providing evidence of a competitively structured community. Secondly, the V-ratio index for the two algorithms also showed a non-random pattern of species association. V-ratio observed values (0.7; p value < 0.05) were significantly smaller than those expected by chance, showing negative covariance between species pairs, being evidence again of a competitively structured community. Finally, for the two algorithms, the C-Score SES was larger than 2 and smaller than − 2 for the V-ratio, providing evidence of non-random community patterns (Table 1). Contingency table analysis indicated that TCM was detected in 162 (59%) of the 275 orchards sampled and was absent in 113 (41%). There was competition between E. banksi and E. orientalis, since E. orientalis was only present in 30 (19%) of the orchards where TCM occurred, while it was absent in the remaining 132 plots (81%). In the orchards where TCM was absent, E. orientalis increased its presence (38; 34%) and decreased its absence (75; 66%) (χ 2 = 7.37, df = 1, p value = 0.0066) (Fig. 2a). Furthermore, a competitive relationship between E. banksi and P. citri was also detected (Fig. 2b). In the same way as in the previous case, the presence of P. citri increased in the orchards where TCM was absent (49; 43%) and decreased with its presence (47; 29%). In parallel, its absence also gained strength with the presence (115; 71%) of TCM and declined with the absence (64; 57%) (χ 2 = 5.42, df = 1, p value = 0.019). On the other hand, no statistically significant evidence of a competitive interaction between E. orientalis and P. citri could be found (Fig. 2c) (χ 2 = 0.0048, df = 1, p value = 0.9443). Probabilistic model of species co-occurrence also allowed us to analyse the co-occurrence between pairs of species in the mite community. Eutetranychus banksi and E. orientalis significantly coexisted (30) in fewer orchards than expected by chance (40.1) showing a strong competitive relationship (p-lt = 0.0034). In the same way, TCM had a negative relationship with P. citri as both also appear significantly together in fewer orchards (47) than expected (56.6) (p-lt = 0.0101). However, between P. citri and E. orientalis it was not possible to detect either a negative (p-lt = 0.4751) or a positive (p-gt = 0.6391) relationship. The effect size of the competition for E. banksi with E. orientalis and P. citri was similar (Table 2).

Species distribution and GLM competition models
The current geographic distribution of the three mite species is shown in Fig. 3. TCM showed the highest mean density values with 0.0068 occurrences/km 2 distributed in three high-density areas with about 0.05 occurrences/km 2 . The most important of these areas is located in the south of the Valencia province. Furthermore, two other smaller areas can also be observed, the first between the provinces of Valencia and Castellón, and the second one between the provinces of Valencia and Alicante (Fig. 3 left). Panonychus citri presented 0.0041 occurrences/km 2 of mean density with two high-density areas with around 0.05 occurrences/km 2 . The most important of them coincides with the northern nucleus of TCM between the provinces of Castellón and Valencia, and the second smaller one is in the north of the province of Castellón (Fig. 3 middle). Finally, E. orientalis with 0.0028 occurrences/km 2 showed the lowest mean density value in two low density areas (0.026 occurrences/km 2 ). The first one located in the north The mean TCM density in orchards without E. orientalis (0.028 ± 9.59e − 04) was significantly higher than in orchards with E. orientalis (0.016 ± 1.55e − 03) (W = 10,470, p value = 1.631e − 09) (Fig. 4a, left). To a lesser extent, P. citri also showed the highest TCM density values in orchards with its absence (0.026 ± 1.08e − 03), while orchards with its presence showed low values (0.022 ± 1.45e − 03) (W = 10,113, p value = 0.01558) (Fig. 4a, right). The mean density values of E. orientalis was significantly higher in orchards not infested with TCM (0.010 ± 6.95e − 04) and lower in infested ones (0.006 ± 5.56e − 04) (W = 11,845, p value = 3.355e − 05) (Fig. 4b, left). There was no significant difference for the E. orientalis density values between orchards with or without P. citri (W = 8040, p value = 0.3804) (Fig. 4b, right). Finally, P. citri mean density value was significantly lower in the presence (0.016 ± 1.01e − 03) of TCM than in its absence (0.022 ± 1.43e − 03) (W = 11,262, p value = 0.001159) (Fig. 4c, left), and it was not significantly affected by the presence of E. orientalis (W = 7294, p value = 0.6534) (Fig. 4c, right). A summary for the GLM analysis is provided in Table 3. The first model for TCM shows a significant steep decline in the presence of E. orientalis with increasing E. banksi density (Fig. 4g) (Table 3a). To a lesser extent, the presence of P. citri was also negatively affected by the density of E. banksi (Fig. 4d) (Table 3b). Complementarily, in E. orientalis and P. citri models, the probability of finding TCM was strongly affected by the density of E. orientalis (Fig. 4e), and in a smaller proportion by P. citri (Fig. 4f) (Table 3c). There was no statistically significant evidence about the influence of E. orientalis density on the presence of P. citri and vice versa.

Multivariate analysis
PCA shows a clear difference between the orchards colonised by each of the mite species. The first axis explains 65.72% of the total variance, while the second axis explains 34.28% (Fig. 5). PERMANOVA test provided statistical support for the visual interpretation of the ordination results, and evidence that samples within groups are more similar than Table 1 Results for the Null model co-occurrence analysis. The analysis was based on the C-Score and V-ratio index for sim 2 and sim 10 algorithms Competition samples from different groups. Differences were highly significant, due only in a small proportion to multivariate dispersion (Table 4a, b). Furthermore, PERMANOVA models including only pairs of mite species also indicated that there were highly significant differences between groups (Table 4c). Across all the sampled locations with at least 1 mite species (241), the initial 17 environmental variables were reduced to 12 using VIF < 10, and FS for environmental RDA did not exclude any variable Fig. 2 Proportion of orchards with occurrence of each species (framed), when another one occurs (grey legend) or is absent (white legend). a E. orientalis versus E. banksi, b P. citri ver-sus E. banksi, c P. citri versus E. orientalis. For each pair, significant differences are denoted with asterisks. Chi square contingency test: *P < 0.05; **P < 0.01 Table 2 Results for the probabilistic model of co-occurrence analysis Signif. Codes: *p < 0.05, **p < 0.01 (Obs-cooccur) = Observed number of orchards in which the species overlap; (Exp-cooccur) = Expected number of orchards in which the species should coincide under random distribution; (p-lt) = Probability that two mite species would co-occur less frequent than the observed number of co-occurrences. It can be interpreted as p value indicating significance level for negative co-occurrence patterns with significance threshold at p < 0.05; (p-gt) = Probability that two mite species would co-occur more frequent than the observed number of co-occurrences. It can be interpreted as p value indicating significance level for positive co-occurrence patterns with significance threshold at p < 0.05. (Effect size) = Indicates the intensity of the relationship, negative (x < 0), positive (x > 0) or random (x = 0) The RDA constrained by the selected environmental variables across all the sampled locations ( Fig. 1 in Supplementary Material 5) indicated that TCM is associated with high temperatures as well as not being negatively affected by rainfall or high relative humidity. On the other hand, P. citri is found in low temperatures, high relative humidity, and areas with high wind speed, and E. orientalis is found in areas with high temperatures, slow wind velocity and low relative humidity. TCM was present in areas with the highest density of citrus crops. A summary of the environmental variables for the sites where each mite species was found is shown in Table 2 in Supplementary Material 2. The results for the RDAs constrained by the selected environmental and spatial variables, and VPA are summarised in Table 5. Selected environmental variables explain 75.4% of the total observed variance across all sampled locations. However, despite this, most of the variance explained by the environment is shared with space (74.9%). If we look at the variance explained purely by the environment (0.5%), it is much smaller than the variance explained purely by space (22.2%). The shared variance is also high for E. banksi, E. orientalis and P. citri sites. Furthermore, the spatial fraction was also more important than the environmental fraction in all three scenarios. However, E. orientalis with 64.6% showed the highest purely spatial fraction followed by E. banksi (46.3%) and finally P. citri (23.8%). On the other hand, P. citri with 1% had the best purely environmental fraction followed by E. banksi (0.9%) and E. orientalis (0.1%) ( Table 5).

Discussion
Agricultural ecosystems have become excellent laboratories for the study of invasion ecology. The increased international mobility of people and goods during the Table 3 Results of the binomial (logit-link) generalized linear models (GLMs) on the probability of presence depending on the density Kernel (KDE) for the three mite species Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0. recent decades has caused the accidental introduction of invasive pests in agriculture to rise, leading to economic losses for farmers and food producers (Pimentel et al. 2000). Furthermore, the ecological impact produced by the invaders may change the resident community, altering the delicate web of relationships among species (Parker et al. 1999). Phytophagous mites have high potential as an invasive species because their small size (from 0.2 to 0.5 mm length), short generation time, high reproductive potential, and their ability to spread by the wind. Despite this, the study of the ecological consequences after invasion and the competitive relationships between spider mite species, the most damaging group to cultivated plants, has only been carried out on a few cases (Reitz and Trumble 2002). Recently, Ferragut et al. (2013a) studied the mite species composition on cultivated and non-cultivated plants in vegetable agrosystems before and after the invasion by the tomato spider mite Tetranychus evansi Baker and Pritchard, demonstrating that the invasive mite had reduced the presence and abundance of other Tetranychus species. However, a different approach was applied in that study and no spatial analysis was carried out to assess interspecific competition. The present work is the first study to conduct a spatial analysis to examine interspecific competition among spider mite species.

Occurrence analysis
We have studied the impact of TCM after its arrival. Our results have confirmed the rapid spread of this species since 2013, becoming the most frequent and abundant spider mite on citrus and leading to a significant reduction in the presence of their relatives P. citri and E. orientalis. These conclusions are supported

Species distribution and GLM competition models
The current geographic distribution of the three spider mites do not coincide geographically, and GLM analysis shows how the probability of finding E. orientalis or P. citri decreases with increasing TCM density and vice versa, evidence of current competitive relationships. Unfortunately, there is no information about the geographic distribution of E. orientalis prior to the invasion of TCM. On the other hand, P. citri spread rapidly to the north and south of the entire Valencian territory soon after the first detection in April 1981, causing higher incidence in the Valencia province (García-Marí et al. 1983). A survey carried out by Abad-Moyano et al. (2009) in citrus crops, before the arrival of TCM, clearly showed P. citri widely distributed throughout the citrus-growing area of the provinces of Castellón and Valencia. According to our results, the arrival of TCM to the province of Valencia has decreased the presence of P. citri, which now has its core distribution further north, in the south of the province of Castellón.

Multivariate analysis
PCA and PERMANOVA indicated that there were segregated geographical distribution patterns between all pairs of mite species. Segregated patterns can be explained by competitive relationships between species. Nevertheless, they also may be caused by abiotic factors such as environmental heterogeneity across a set of sites (Kraft et al. 2015), acting as filters that hinder the establishment of the species (Pearson et al. 2018). Despite the economic importance of spider mites, there are no studies on the influence of environmental factors on the structure of their communities, and climatic factors such as temperature, relative humidity (RH), and precipitation have been considered the most influential in the development of their populations (van de Vrie et al. 1972;Vacante 2010). Furthermore, wind velocity or crop density may also be underlying non-random distribution patterns, since spider mites are considered passive dispersers, and wind dispersal has been proposed as one means of long-distance dissemination, especially when spider mite colonies are overcrowded, or the host plant is no longer suitable (Boyle 1956;van de Vrie et al. 1972). VPA allows partitioning of the variance of species abundance community data into variance that can be explained by environmental factors or by spatial variables. In practice, in a VPA the spatial component represents the community variance that cannot be explained by the environmental factors considered. In our study we have selected two sets of variables, a first set of environmental variables related to the temperature, relative humidity, precipitation, wind speed and crop density; and a second set of spatial variables that tries to integrate the spatial variance of the set of sites. RDA constrained by the environmental variables across all the sampled sites suggests that the spider mite communities were significantly affected by these environmental factors, explaining 75.4% of the community variance. However, when we separate the variance due to the environment, as well as the shared variance using VPA, the purely environmental contribution to the variance of mite communities (0.5%) is much lower than the spatial one (22.2%) (non-environmental). Furthermore, we found a large proportion of the explained variance shared by the environmental and spatial variables (74.9%). This situation is common, since environmental factors often vary through space, especially on a broad scale. In this scenario, we should be careful when assigning cause-effect relationships between species and space or with environmental conditions (Borcard et al. 1992). Furthermore, purely spatial fraction also could be explained by unmeasured environmental variables (Borcard and Legendre 1994). Despite these two observations, our findings indicate the existence of spatial patterns of geographical distribution of these species which do not seem to be governed mainly by the environmental variables selected. These results agree with Ferragut et al. (2013b) who point out that the climatic conditions under which citrus crops are grown throughout the Valencian region are suitable for the establishment of E. banksi, E. orientalis and P. citri populations all year-round. In addition to the environmental conditions, historical factors and interspecific competition can act to spatially structure the community in VPA (Borcard and Legendre 1994). In our work we have demonstrated the existence of competition among the three spider mite species using different approaches, and most of the classical publications point to competitive interactions as the main important factor in species segregation (Diamond 1975;Connor and Simberloff 1979) determining the community assembly (Kunstler et al. 2012). Furthermore, processes such as colonization/extinction history may also be underlying nonrandom distribution patterns (Gilpin and Diamond 1982). The current geographical distribution for the three mite species is not coincidental and evidences their competitive relationships. Moreover, the density cores of the recently introduced TCM and E. orientalis are related to the place of their first detection, the south of the province of Valencia, and the south of the province of Alicante, respectively. For all these reasons, we believe that competitive interactions, and the colonization history, are the most important factors conditioning the spider mite community structure across all the sampled sites, contributing to the spatial component in the VPA. In our VPA study for each species, we found that (as in the analysis for all the sampled localities), the pure spatial component is much greater than the pure environmental component. If we focus on the analysis of the purely environmental fraction, P. citri followed by TCM have the highest values, while E. orientalis has a very low value. However, if we pay attention to the purely spatial component, E. orientalis has the highest value, followed by E. banksi and finally P. citri. We believe that these results, can be explained by historical colonisation factors and competitive interactions between species. According to the historical colonisation, P. citri has been present in the region for about 40 years; therefore, there is currently no spatial evidence of its colonisation process (historical factors) in its spatial distribution, and its pure environmental component is large, reflecting its adaptive process. Despite this, its current competition with TCM is reflected in its pure spatial component. Eutetranychus banksi, although recently introduced, has spread rapidly throughout the citrus area displacing the other two mite species and, accordingly, its spatial distribution is less affected by competition and historical colonisation processes (spatial component) compared to E. orientalis. However, its rapid adaptation process is reflected in its broad environmental component. Finally, E. orientalis has a very small pure environmental component, as it has a very large pure spatial component due to intense competition with TCM and the historical colonisation processes.

Competitive mechanisms
In our opinion, there are different non-exclusive hypotheses that can explain the impact of TCM on the whole community. One of the best-known hypotheses to explain the spread of invasive species is the Enemy Release Hypothesis. Following this approach, invasive species are less impacted by natural enemies (e.g., predators, parasites, and pathogens) than native species, because, in the new geographic location, invasive species lack the enemies that keep their growth under control in their native environment (Jeffries and Lawton 1984). There is little information available on field efficacy of predators on TCM in its native geographical area. The predatory mites of the family Phytoseiidae are the most important biological control agents used against spider mites, and essential elements in the pest management programs (McMurtry 1982). In Florida, the phytoseiids Euseius mesembrinus (Dean) (Abou-Setta and Childers 1989), Iphiseiodes quadripilis (Banks) (Villanueva and Childers 2007) and Galendromus helveolus (Chant) (Caceres and Childers 1991) are common species on citrus and were tested in the laboratory with TCM as a food source, being a suitable prey. Furthermore, Landeros et al. (2004) observed a synchronization of the population peaks of TCM with their predator E. mesembrinus, as well as a correlation in the abundance of both species on Mexican citrus crops. Nevertheless, these natural enemies were not introduced with TCM in the Spanish invaded areas. Euseius stipulatus (Athias-Henriot) is the predominant native phytoseiid on citrus in Spain and is a useful biocontrol agent of P. citri (Ferragut et al. 1987(Ferragut et al. , 1988. According to McMurtry (1985), predators that attack P. citri may also be effective against Eutetranychus spp. because those spider mites have a similar colonization pattern on the leaves and produce small amounts of webbing. This predator was studied under laboratory conditions to test its predatory potential on E. banksi. The predator was able to complete its life cycle when feeding on TCM, but mortality was high, the reproductive parameters very low and almost no eggs were produced, indicating that this prey is not a suitable food source (data not published). According to the results obtained in laboratory tests, we believe that E. stipulatus, in the presence of P. citri and E. banksi as a food source, will prefer the former. Therefore, E. stipulatus, by selecting P. citri over E. banksi, would be contributing to tip the balance in favour of TCM in the interspecific competition between both phytophagous mites. Despite the spread of TCM, E. stipulatus populations do not seem to become reduced, since this phytoseiid is highly generalist and can complete its development feeding on other mite species, small insects, and even pollen or fungal spores (Ferragut et al. 1987;McMurtry and Croft 1997).
In addition to the ability of TCM to escape the attack by native predators, the biology of this species can turn it into a superior competitor of P. citri at high temperatures. The developmental time of TCM at 25 °C is 13.1 days , while at the same temperature P. citri has a shorter development time of 10.9 days (Tanaka and Inoue 1970). On the contrary, at higher temperatures the development time of TCM is considerably reduced, being of 11.6 days at 28 °C and 61% RH. Greater differences are reflected in the fecundity values. The highest fecundity for TCM was reported at 28 °C with 37.08 eggs/female, with 8.84 being the maximum fecundity rate (eggs/female/day) at 30 °C. This behaviour is the opposite of that of P. citri, which produces 40.4 eggs/female and 6 eggs/female/day at 24 °C (Yasuda 1982). The temperature range of 28-30 °C seems to be optimal for TCM, while the range of 24-25 °C is optimal for P. citri. In the Valencian region in late spring-early summer, with average daily temperatures between 24-26 °C, low populations of P. citri can be detected which are controlled by E. stipulatus. During July-September, mean daily temperatures reach values of 27-30 °C, providing conditions more favourable for the development of TCM. Low populations of this mite are found in July, increasing rapidly to produce population peaks at the end of September and beginning of October (data not published). According to Ferragut et al. (1988), P. citri has a population peak in late summer-early fall. However, this peak is currently not observed in the field (data not published). We believe that the demise of the population peak of P. citri is related to the lack of effective natural enemies against TCM, as well as the indirect competition for resource exploitation between both species. TCM seems not to reject feeding on leaves slightly damaged by previous consumption of low P. citri populations, while high TCM populations completely cover the habitable space on the leaves and leave them without nutrients available for P. citri. The Resource Depletion hypothesis (space and nutrient) has been cited as the leading mechanism of interspecific competition between the spider mites Tetranychus urticae Koch and Eotetranychus carpini borealis (Ewing) on red raspberry (Bounfour and Tanigoshi 2001), while the ability of Tetranychus telarius (L.) to feed on leaves previously consumed by Panonychus ulmi (Koch) gives them an advantage in a competitive situation on peach and apple crops (Foott 1963).
In addition to these two hypotheses, some studies have shown how indirect competition can exist between two herbivores through changes in the quality of the shared plant, representing interactions between competitors mediated by the host plant (Karban and Myers 1989;Tallamy and Raupp 1991 (Sarmento et al. 2011). Therefore, the induction of plant defences and their tolerance by the inducing species determine the outcome of their competition with less tolerant species. Although there are no studies on the ability of TCM to manipulate citrus defences in its favour to affect potential competitors, we cannot rule out that this mechanism might play a role in the competitive displacement carried out by E. banksi.
Our results show, for the first time, the existence of competition between two invasive mite species of the genus Eutetranychus originating from different parts of the world. In this case, we cannot consider the enemy release hypothesis as a mechanism that regulates the competition because E. orientalis, like E. banksi, lacks effective predators in the Spanish citrus area. However, our results point to the existence of a competitive relationship between them that is even more intense than that recorded between E. banksi and P. citri. Both species are phylogenetically closer, so their ecological niche may be more similar, and therefore the indirect competition through resource exploitation between both species would be more intense. Previous laboratory studies demonstrated that E. orientalis performs better at high temperatures around 30 °C, in the same way as TCM; however, its biological parameters, especially those related to fecundity, are lower than those of TCM at the same environmental conditions. Eutetranychus orientalis showed a fecundity of 14.56 eggs/female at 25 °C and 16.33 eggs/female at 30 °C (Imani and Shishehbor 2009), while TCM can produce twice as many eggs as E. orientalis (29.96 eggs/female at 25 °C and 32.08 eggs/female at 30 °C) . Fecundity Differences between both species would be able to lead the competitive displacement of E. orientalis from the citrus area (Reitz and Trumble 2002).
As both Eutetranychus are closely related species, Reproductive Interference is another hypothesis that may explain the intense competition. Reproductive interference is an interspecific interaction during the process of mate acquisition caused by incomplete species recognition resulting in the lack of courtship and mating discrimination, which adversely affects the fitness of at least one of the species involved (Gröning and Hochkirch 2008). Since the native ranges of both Eutetranychus do not overlap, they likely do not have recognition mechanisms that act as a barrier to avoid interspecific mating. Geographic displacement driven by this mechanism has been demonstrated in spider mites on apple trees in Japan, where Panonychus mori Yokoyama predominates in the north and P. citri in the south. Interspecific mating occurs without fertilization and, because of their reproduction mode by arrhenotoky, females involved in interspecific copulation do not produce females in their offspring, resulting in a decrease in their biological fitness (Fujimoto et al. 1996). However, the degree of reproductive interference is greater in P. mori than in P. citri and, consequently, P. citri geographically has displaced P. mori, in the southern area (Takafuji et al. 1997). Following this hypothesis, the intense competition detected by us between E. banksi and E. orientalis could suggest the existence of reproductive interference between them. It is possible that males of E. orientalis do not discriminate between conspecific and heterospecific females as much as males of TCM, or alternatively, females of TCM discriminate against interspecific mating more than females of E. orientalis. In either scenario, E. orientalis females may be involved in interspecific mating more frequently than TCM females, being more affected by reproductive interference. Interspecific copulations could be totally ineffective, but not affect the successive fertilization with males of the same species (Ozawa and Takafuji 1987). However, in the worst case, the eggs produced are non-viable, no females are produced in the offspring, or hybrid and infertile females are generated, affecting in all cases subsequent intraspecific mating (Boudreaux 1963;Smith 1975).
Niche partitioning and future scenarios Interspecific competition may cause competitors to evolve to reduce the impact of that competition (Sakai et al. 2001). We can consider future scenarios at the plant-leaf level and at plant-species level. One way to reduce the impact of competition and make possible the coexistence between competitive species should involve a niche partitioning process (Albrecht and Gotelli 2001). In a broad scale spatial niche partitioning, P. citri could be displaced to more temperate citrus-growing areas, where it can compete better, specializing in those temperate ranges, while TCM would occupy the warmer zones. However, this option seems less plausible in the case of E. orientalis, as both Eutetranychus prefer hot temperatures. On the other hand, on a fine scale, the movement of inferior competitors E. orientalis and P. citri to the lower face of the citrus leaves avoiding the presence of TCM, who prefers to live on the upper side, could allow the coexistence of these two species on the same leaves with their superior competitor. However, this hypothesis seems unlikely, at least in clementine citrus and lemons, since the lower side is already occupied by the two-spotted spider mite T. urticae, which produces large amounts of silk covering the leaf surface and probably rendering this habitat unfavourable (Morimoto et al. 2006).
A food niche partitioning event was observed after the invasion of T. evansi in the Valencian region, as the native spider mite species T. urticae and T. turkestani Koch moved to occupy non-crop plants more frequently, while T. evansi established mainly on crop plants (Ferragut et al. 2013a). Thus, another option should be a change in the host plant range by the inferior competitors in the invaded area. The three mite species are highly polyphagous, reported on more than one hundred plant species worldwide (Migeon and Dorkeld 2022) and the observed interactions may be expected on other hosts. We have demonstrated the existence of competition on citrus plants, but there is no information on the population behaviour of these mites on other plants.
Although food and space are the major mechanisms for niche partitioning, time can also provide a scenario for the niche partitioning on an annual temporal scale (Loureau 1989). This scenario is not very plausible for both Eutetranychus species that show optimal biological parameters related to high temperatures. However, P. citri could present a spring population dynamic on warm temperatures, while TCM populations could develop during the hot summer. Nevertheless, this scenario is improbable, since E. stipulatus, the main predator of P. citri, has its population peak in spring (Ferragut et al. 1988). Competition among E. banksi and P. citri in Florida citrus has been reported before by Childers (1992), where TCM has become the dominant species in the orchards since 1955, displacing the P. citri. Therefore, the establishment of E. banksi and the resulting competitive displacement and extinction of P. citri from entire or most of the citrus growing areas in the Valencian region seems the most likely hypothesis. On the other hand, this is the first time that E. orientalis and E. banksi have coexisted in citrus crops worldwide, but the phylogenetic proximity, and the consequent niche similarity and intense competition between both Eutetranychus species, makes their coexistence even less plausible. Whatever the final outcome, under a total competitive displacement scenario, TCM will be more exposed to intraspecific competition, a powerful force in the case of agricultural pests which, by definition, produce intense population outbreaks in short periods of time.