Global patterns and drivers of herbivorous eriophyoid mite species diversity

Environmental drivers and host richness play key roles in affecting herbivore diversity. However, the relative effects of these factors and their effects on lineages characterized by high host specificity are not well known. In this study, we explored the extent to which contemporary climate, Quaternary climate change, habitat heterogeneity and host plants determine the species richness and endemism patterns of herbivorous eriophyoid mites.


| INTRODUC TI ON
Plants and herbivores comprise the highest macroscopic terrestrial diversity (Futuyma & Agrawal, 2009;Pacala & Crawley, 1992), but they are unevenly distributed at the geographical scale (Myers et al., 2000). This uneven distribution and its potential drivers have fascinated ecologists and biogeographers for over two centuries (Araújo et al., 2008;Humboldt & Bonpland, 1807;Martins et al., 2021;Mittelbach et al., 2007;Rahbek et al., 2019). Several hypotheses have been proposed in the last two decades to explain these biodiversity patterns (Table 1). Clearly, the current distribution of species diversity is shaped by complex interactions, including biotic and abiotic factors (Hawkins et al., 2003;Svenning & Skov, 2005). However, the specific impacts of biotic and abiotic factors in shaping the current phytophagous mite species diversity are not well understood.
Eriophyoid mites are of particular interest because they comprise the most species-rich superfamily (Eriophyoidea) in the Acari (ca. 4400 species), covering more than 60% of phytophagous mites (Zhang, 2011). They exhibit a worldwide distribution and be divided into three families: Phytoptidae (ca. 160 species), Eriophyidae (ca. 3790 species) and Diptilomiopidae (ca. 450 species) (Amrine Jr. et al., 2003;Zhang, 2011). All eriophyoid mite species are strictly phytophagous, resulting in intricate relationships with host plants, such as by making galls or blisters (known as gall mites), curling leaves, transmitting plant viruses or living as vagrants on leaf surfaces; therefore, some of these species are considered as pests in agriculture and forestry and can cause massive economic losses (de Lillo et al., 2018). Besides, eriophyoid mites can be easily distinguished from other mites by their small size (~200 μm in length) and by having only two pairs of legs (known as four-legged mites) (Amrine Jr. et al., 2003).
Numerous studies on eriophyoid mites have been undertaken over the past 150 years (Amrine Jr. et al., 2003), but the mechanisms underlying their diversity remain poorly understood.
Eriophyoid mites have high host specificity, that is, over 80% of species are monophagous . Therefore, it is possible that host plants determine their global richness and drive their endemic species richness, which might in line with the 'resource dependent hypothesis' (Du et al., 2020;Futuyma & Agrawal, 2009). Moreover, they are unevenly distributed on a global scale (Amrine Jr. & Stasny, 1994). Previous studies pointed that they are distributed mainly in temperate regions and occur less in tropical and subtropical regions (de Lillo & Skoracka, 2010).
Because eriophyoid mites are strictly phytophagous, the occurrence of high host plant richness could provide more resources for herbivores (Futuyma & Agrawal, 2009). Therefore, high eriophyoid mite richness would be expected in tropical regions.
However, no strict analyses based on comprehensive datasets have been performed to examine their global species richness patterns. Alongside factors of host plant, other factors that may contribute to their diversity include historical climate fluctuations (e.g., Quaternary climate change), the contemporary climate (e.g., environmental energy, water availability and climate seasonality) and habitat heterogeneity (i.e., topographic heterogeneity, percentage of global tree coverage).
Global species diversity reflects the species richness and endemic species richness (Du et al., 2020;Myers et al., 2000;Phillips et al., 2019;Possingham & Wilson, 2005). In this study, we compiled a dataset covering 97% of extant eriophyoid mite species (ca. 4278 species) and 22,973 occurrence sites to (1) uncover the species richness and areas of endemism (AoE) patterns of these mites across six zoogeographical regions (Figure 1a; Wallace, 1876), (2) Figure S1). The occurrence sites were determined to the best possible extent, thus forming a large dataset including almost all extant eriophyoid mite records and spanning a wide spectrum of vegetation types.
Records of eriophyoid mite occurrence and their host plant are abundant across the Northern Hemisphere, but relatively sparse in the tropics. To avoid the impact of sampling bias on statistical analysis, we plotted rarefaction curves based on the total number of eriophyoid mite species and host plant species, respectively, per 300 × 300 km grid using the iNEXT function in R (Hsieh et al., 2016).
We counted the number of individuals for each species per grid and then used rarefaction/extrapolation (92.5% sample coverage) to extrapolate an unbiased richness value per grid. This unbiased richness value was used for global species richness, endemic species richness and host plant richness analyses. The iNEXT package computes and plots the sample size and coverage-based rarefaction and extrapolation of Hill numbers, which can be rarefied to a smaller sample size or extrapolated to a larger sample size (Chao & Jost, 2012;Gotelli & Colwell, 2001). This approach can standardize the sampling size and can be used to predict the potential species richness (Colwell et al., 2012).

| Eriophyoid mite species richness and endemism
To identify the global eriophyoid mite species richness, we divided the global map into 300 × 300 km grids and mapped the unbiased species richness based on Lambert cylindrical equal-area projection using ArcGIS 10.7 (ESRI, Inc.). The total number of eriophyoid mites in each grid was counted and ranked. The grades were shown TA B L E 1 Predictions and hypotheses for herbivorous species diversity influenced by environmental factors, latitude gradients, host plants and historical climate stability

Graphical prediction Hypothesis
'Species-energy hypothesis' suggests that the available energy in different regions determines species richness (Hawkins et al., 2003;Wright, 1983) 'Biotic interactions hypothesis' predicts that the strength of biotic interactions (herbivore pressure) declines with latitude, leading to a greater diversity of plants and herbivores in low latitude regions (Baskett & Schemske, 2018;Schemske et al., 2009) 'Resource dependent hypothesis' supports that the diversity of herbivore and host plants have intricate relationships, in which herbivores have high host specificity, infesting one or a few related host plants ; thus, host plant resources might constrain specialist herbivore diversity (Du et al., 2020;Futuyma & Agrawal, 2009) 'Heterogeneity hypothesis' holds that suitable habitats can provide more niche opportunities and diverse topography for species colonization, speciation and diversification (Stein et al., 2014(Stein et al., , 2015 'Climate stability hypothesis' concludes that long-term climate fluctuations have different effects on species diversity with different dispersal capacities, while stable climate increases net speciation by reducing extinction rate (Fine, 2015;Mckenna & Farrell, 2006;Sandel et al., 2011) in different colours on the global map. The top 7.5% of richest cells were identified as hotspots.
Areas of endemism (AoEs) were analysed using the heuristic algorithm of NDM/VNDM (version 3.1) (Prado et al., 2014;Szumik & Goloboff, 2004). NDM/VNDM determines the AoEs based on the optimal criteria by assigning an endemicity value to each grid regardless of how that area was found or hypothesized and considering the spatial position of taxa as an integral part of the analysis. We mapped the AoEs to latitude-longitude grid sizes of 2° (220 × 220 km), 3° (330 × 330 km) and 4° (440 × 440 km) using Lambert cylindrical equal-area projection. It has been suggested that the use of three different grid sizes improves the reliability of inferred endemism patterns (Aagesen et al., 2009;Escalante et al., 2010). Areas with a high degree of overlap for different grid sizes reflect a high degree of endemicity (Prado et al., 2014). We input latitude and longitude data for each species to compute endemicity. Scores over 2.0 were saved for each grid which indicate that two species coincide in their distributional areas. Then, we examined subsets with 50% unique species as overlapping and searched for 50 times. Finally, we selected species with a minimum score of 0.4 and calculated the consensus areas using a cut-off with 100% similarity.

| Environmental variables and other factors
To test the relative effects of abiotic variables on global eriophyoid mite species diversity patterns, we selected 24 environmental variables as potential drivers and grouped them into three themes, namely contemporary climate, Quaternary climate change and habitat heterogeneity (Table S2). Environmental energy, water availability and contemporary climate seasonality were chosen to represent contemporary climate. The Quaternary climate change data included climate of the last glacial maximum (LGM; ca. 21,000 years ago) and anomalies representing long-term climate change between the LGM and current climate (Liu et al., 2017;Sandel et al., 2011).
The elevation range (topographic heterogeneity, TOPO) and vegetation (percentage of global tree coverage, PTC) were selected as habitat heterogeneity factors that influence species richness at broad spatial scales.
The contemporary climatic variables were obtained from the WorldClim dataset (http://www.world clim.com/version2) at a spatial resolution of 30 arc seconds (~1 km 2 ), from 1970 to 2000 (Stevens, 1989). Quaternary climate change (annual mean temperature during the last glacial maximum, LGM TEMP) and annual mean precipitation during the LGM (LGM PREC) were also downloaded io/) at a spatial resolution of 30 arc seconds. The key ecological variables were extracted and resampled for subsequent analyses based on 300 km × 300 km grid cells in ArcGIS 10.7.

| Statistical analyses
Assessing the impacts of environmental factors on species richness is controversial because spatial autocorrelations may exist between the occurrence sites and nearby spatial locations (Dormann et al., 2013;Hawkins, 2012). This can lead to a misinterpretation between the explanatory variables and species richness (Kubota et al., 2015). Therefore, it is crucial to avoid multicollinearity among climate variables and their spatial autocorrelations (SACs) in the context of assessing the effects of environmental variables on species richness. We used multiple approaches to avoid these shortcomings. All variables were standardized (using z-scores) to eliminate dimensional effects between ecological factors. We then removed the highly correlated variables (|r| > 0.7) (Dormann et al., 2013).
Finally, we tested the collinearity between variables using a variance inflation factor (VIF) with a threshold of less than three (Zuur et al., 2010). Six of the 24 environmental factors were retained for simultaneous autoregressive (SAR) and ordinary least squares (OLS) regression analysis (Tables S3 and S4). We further added host plant richness as one additional factor for SAR and OLS analyses.
The SAR model and common OLS regression method have been widely used in multiple regression analyses, which accounts for the association between species richness and contemporary climate, Quaternary climate change and habitat heterogeneity (Araújo et al., 2008;Du et al., 2020;Kissling & Carl, 2007). SAR models are more adequate than other spatial regression approaches (Beale et al., 2010) because they can be used to integrate spatial dependence into error terms. We calculated the best optimal combined model from 367 models according to Moran's I value, together with the Akaike information criterion (AIC) and maximum model fit (R 2 ).
We used a random forest model to examine the influence of different sets of drivers on species diversity using R package 'randomforest' 4.6-14 (Breiman, 2001). The random forest model is less sensitive to spatial autocorrelation (Sinha et al., 2019); meanwhile, the nonlinear relationship between species diversity and various environment variables can be examined (Li et al., 2021).

| Global diversity patterns of eriophyoid mites
Our

| Factors affecting eriophyoid mite species diversity
Our ordinary least squares (OLS) and simultaneous autoregressive   (Tables S5-S7). The random forest analyses showed that host plant richness is the best explanatory variable for eriophyoid mite richness (R 2weighted importance = 54.132%), followed by Quaternary climate change (LGM TEMP, R 2 -weighted importance = 5.948%; Pano, R 2weighted importance = 3.635%, Figure 4a). The other factors had relatively little effects (Figure 4a). Host plant richness also had the best explanatory variable when accounting for families ( Figure S3).  Table S4). Random forest results also showed that host plant richness was the strongest explanatory variable for endemic species richness (R 2 weighted importance = 49.120%, Figure 4b), followed by Quaternary climate change (LGM TEMP, R 2 weighted importance = 9.912%; Pano, R 2 weighted importance = 4.145%,   LGM PREC, annual mean precipitation during last glacial maximum;

| Biotic and abiotic effects on eriophyoid mite
LGM TEMP, annual mean temperature during last glacial maximum; Pano, absolute anomaly in annual mean precipitation between last glacial maximum and current climate, PTC, percentage of global tree coverage. See Table S2 for definitions of abbreviations.  Table S2 for definitions of abbreviations. ***p < 0.001, **p < 0.01, *p < 0.05. LGM PREC, annual mean precipitation during last glacial maximum; LGM TEMP, annual mean temperature during last glacial maximum; Pano, absolute anomaly in annual mean precipitation between last glacial maximum and current climate; PTC, percentage of global tree coverage.

Contemporary climate
Quaternary climate change Habitat heterogeneity -0.148*

| DISCUSS ION
Our results reveal a pattern in which eriophyoid mite species richness and endemic species richness are higher in temperate regions, opposite to the latitudinal diversity gradient proposed by 'biotic interactions hypothesis' (Baskett & Schemske, 2018;Schemske et al., 2009). We find a driver of climate factors and heterogeneity factors interfered diversity via host plants, in agreement with the 'climate stability hypothesis' (Fine, 2015;Mckenna & Farrell, 2006;Sandel et al., 2011), 'resource-dependent hypothesis' (Du et al., 2020;Futuyma & Agrawal, 2009) and 'heterogeneity hypothesis' (Stein et al., 2014(Stein et al., , 2015. Contemporary climate (annual mean precipitation) was negatively correlated with eriophyoid mite species richness, which is in contrast with the 'species-energy hypothesis' (Hawkins et al., 2003;Wright, 1983). While little is known about the global patterns of herbivorous mite species diversity, our result that eriophyoid mite species richness peaks in the temperate regions is congruent with findings in specialist drepanosiphine aphids (Du et al., 2020), gall-inducing insects (Price et al., 1998) and earthworms (Phillips et al., 2019).

| Species richness and endemism of eriophyoid mites peak in temperate regions
Nevertheless, this uneven distribution of eriophyoid mite richness covered seven of 25 biodiversity hotspots (Myers et al., 2000), including the California Floristic Province, Mediterranean Basin, Indo-Burma, South-Central China, New Zealand, the Philippines and Brazil's Atlantic Forest. Our results further showed that the species richness centre differed among the Phytoptidae, Eriophyidae and Diptilomiopidae (Figure 1), which may be caused by their host usage preference. For instance, roughly half (84/173) of phytoptid mites infest gymnosperms (Table S1) that are predominantly distributed in the high elevations and high latitudes in the northern hemisphere (Bond, 1989). In contrast, only a few eriophyid mites (109/3584) and few diptilomiopid mites (5/521) infest gymnosperms (Table S1) islands, which act as refugia and centres of speciation (cradle) and regions of maintenance of biodiversity (museum) (Jetz et al., 2004;Orme et al., 2005;Terborgh, 1992). Our AoEs' results further suggest that Europe is the cradle for phytoptid mites (Figure 3b), with multiple cradles for eriophyid mites (Figure 3c), and Asia as well as Europe and South-eastern Brazil are cradles for diptilomiopid mites ( Figure 3d). There is a general agreement that refugia were shaped during tectonic movements and climatic fluctuations in the Last Glacial Maximum of the Quaternary and may play a variety of roles in driving the diversity of plants and invertebrates (Bennett et al., 1991;Ramirez-Barahona & Eguiarte, 2013;Stewart et al., 2010), therefore holding AoEs. There is ample evidence to support the existence of multiple refugia in Southern Europe (Bennett et al., 1991;Taberlet & Cheddadi, 2002;Tzedakis et al., 2002), Eastern and Southeast Asia (López-Pujol et al., 2011), New Zealand (Gardner et al., 2004;Marske et al., 2009) andSouth-eastern Brazil (Ramirez-Barahona &Eguiarte, 2013). Our OLS, SAR and random forest analyses further showed that Quaternary climate change (LGM TEMP) significantly affected species richness and endemism (Tables S3 and S4, Figure 4).
Under these contexts, it is likely that host plants expanded their territory from refugia by long-term adaptive responses to climate change (Comes & Kadereit, 1998;Davis & Shaw, 2001). This would have been accompanied by eriophyoid mites as well as pervasive host-switching (Hoberg & Brooks, 2008), thus shaping their extant distribution.

| Climate factors play key roles in the diversity of eriophyoid mites via host plants
The reasons behind what is driving uneven geographical distribution of species have been of great interest in ecology for many decades (Araújo et al., 2008;Kubota et al., 2015). It is widely assumed that the uneven distribution is dramatically affected by long-term and short-term climatic variations (Hewitt, 1996;Liu et al., 2017;Sandel et al., 2011). Our OLS and SAR analyses indicated that eriophyoid mite richness is significantly affected by the LGM TEMP (p < 0.05, Table S3) and Pano (p < 0.001, Table S3), while endemic species richness is significantly affected by the LGM TEMP (p < 0.001, Table S4). Our random forest analyses showed that the LGM TEMP and Pano have significant associations with species richness and endemic species richness (Figure 4). These results agree with the 'climate stability hypothesis' (Fine, 2015;Mckenna & Farrell, 2006;Sandel et al., 2011), implying that historical variables influenced species survival and radiation (Araújo et al., 2008;Kubota et al., 2015). However, our OLS and SAR analyses showed that annual mean precipitation has a significant negative effect on species richness (Table S3), which is inconsistent with the 'species-energy hypothesis' (Hawkins et al., 2003;Wright, 1983).
High annual mean precipitation should impact not only the survival rate but also short-and long-distance dispersal of eriophyoid mites, thus leading to a negative effect on species richness. The importance of climatic variations in shaping the diversity and distribution patterns across global scale has been consistently identified in plants Kubota et al., 2015), insects (Du et al., 2020;Dunn et al., 2009), earthworms (Phillips et al., 2019) and nematodes . In addition, our OLS and SAR analyses showed that habitat heterogeneity, as shown by percentage of global tree coverage (PTC), has a significant effect on species richness (p < 0.01, Table S3), while elevation has a significant effect on endemic species richness (p < 0.05, Table S4).
These results are consistent with the 'heterogeneity hypothesis' (Stein et al., 2014(Stein et al., , 2015. These findings imply that the contemporary species richness of eriophyoid mites was derived by a complex interplay of environmental factors, including the contemporary climate, Quaternary climate change and habitat heterogeneity. Our OLS, SAR, random forest and SEMs results have consistently shown that host plant richness has a predominately direct impact on the eriophyoid mite richness (Table S3; Figures 4a and 5a) and endemic species richness (Table S4; Figures 4b and 5b), indicating that hosts determine the herbivore diversity. This is consistent with previous studies and predictions of plant-herbivore interactions (Crutsinger et al., 2006;Du et al., 2020) and 'resource dependent hypothesis' (Du et al., 2020;Futuyma & Agrawal, 2009). Eriophyoid mites have a high host specificity  which may be the result of coevolution with host plants (Hoberg & Brooks, 2008). Because most eriophyoid mites do not appear to directly damage to their host plants, we therefore propose that eriophyoid mites have adaptively evolved with their host plants over long-term co-evolution. Co-evolution would presumably reduce the need for plant defence responses. This is supported by the finding that the genome of tomato russet mite has a greatly reduced number of environmental response genes associated with chemoreception and detoxification (Greenhalgh et al., 2020).
Taken together, we suggest that abiotic factors indirectly influenced the richness of eriophyoid mites, via biotic factors-host plants. We therefore highlight complex interactions including biotic and abiotic factors in shaping the extant species richness of specialist herbivores (Hawkins et al., 2003;Svenning & Skov, 2005).

| Limitations and implications
We compiled a dataset comprising almost all records of eriophyoid mites including 22,973 occurrence sites from 102 countries and 4278 species (Figure 1) and identified the key environmental factors (Table S3) that drive this distribution. However, records from some recognized biodiversity hotspots (e.g., Tropical Andes, Central Chile, Sundaland, etc.; Myers et al., 2000) were limited, even with more than 150 years of taxonomic research activities (Amrine Jr. et al., 2003;Amrine Jr. & Stasny, 1994). This might be caused by uneven field surveys or a lack of experts. Future field surveys and faunistic activities are needed to fully understand the global biodiversity of eriophyoid mites. Nevertheless, our research laid a necessary groundwork for future predictions of eriophyoid mite species diversity and underlined some promising strategies for investigating herbivores biodiversity.

| CON CLUS IONS
We synthesized a global scale macroecological dataset and performed multiple regression analyses to decipher herbivore-hostsenvironment relationships. Contemporary climate (annual mean precipitation), Quaternary climate (Pano and LGM TEMP), habitat heterogeneity (PTC) and host plant richness are the main drivers in shaping eriophyoid mite spatial occurrence, while Quaternary climate (LGM TEMP), elevation and host plant richness play key roles in determining endemism distribution. We further demonstrate that the direct effects of climate drivers and environmental drivers are higher on host plant species than on eriophyoid mites, but there is a strong direct positive correlation linking host plants to eriophyoid mites. The current eriophyoid mite richness is shaped by complex interactions, including biotic and abiotic factors. We provide new evidence regarding the global pattern of herbivorous mite diversity.

ACK N OWLED G EM ENTS
This research was funded by the National Natural Science Foundation of China (32170466, 32161143014). The authors thank Liang-Fei Yao and Yue Hu (Department of Entomology, Nanjing Agricultural University, China) for assisting in the data collection on eriophyoid mites. They also thank anonymous reviewers and editors for their thoughtful comments on the manuscript. No permits were required for this study.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used for this study are available from the Dryad Digital Repository (https://doi.org/10.5061/dryad.1ns1r n8v2).