Ecological niche complexity of invasive and native cryptic species of the Bemisia tabaci species complex in China

Bemisia tabaci is an important pest affecting agricultural and horticultural crops worldwide and comprises a complex of cryptic species. In China, the introduction of the two invasive cryptic species, Middle East-Asia Minor 1 (MEAM1) and Mediterranean (MED), has considerably affected the ecological niche of the native cryptic species. Based on occurrence records obtained through field surveys and high-resolution environmental data, using maximum entropy modelling, we established ecological niche models to predict the distribution of invasive and native cryptic species of B. tabaci in China and identified the differences in ecological niches. The results showed that the distribution range and niche breadth of the invasive cryptic species exceed that of the native cryptic species in the order of MED > MEAM1 > China1 > Asia1. There are different degrees of niche overlap and range overlap between cryptic species. Moreover, the important environmental variables affecting their distribution were different, as well as their response and adaptation to most environmental variables. Our results suggest that the B. tabaci species complex occupies a complex ecological niche in China. The findings improve our understanding of the ecological characteristics of B. tabaci species complex, which will be useful in the development of prevention and control strategies for this pest in China.


Introduction
Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae) is one of the most threatening pests in many crops worldwide (De ). This species directly or indirectly affects more than 600 plant species by feeding on phloem, excreting honeydew, and transmitting plant viruses, resulting in significant agricultural losses (Simmons et al. 2008;Jones 2003;Navas-Castillo et al. 2011;Oliveira et al. 2001). Recent genomic studies have shown that B. tabaci is a unique species complex composed of more than 42 morphologically indistinguishable cryptic species (Kanakala and Ghanim 2019;De Barro et al. 2011). Among the cryptic species, the Middle East-Asia Minor 1 (MEAM1, formerly referred to as biotype B) and Mediterranean (MED, formerly referred to as biotype Q) are the two most widely distributed and destructive invasive species (Perring 2001;Tay et al. 2012).
B. tabaci was identified in China in the 1940s, but was not considered a serious pest until the introduction of the two invasive cryptic species. MEAM1 was first identified in China in the mid to late 1990s and has expanded rapidly since then, replacing native cryptic species and Trialeurodes vaporariorum (Hemiptera: Aleyrodidae) in several regions in China (Luo et al. 2002;Li et al. 2016). MED, another invasive cryptic species, was first identified in 2003, and has displaced MEAM1 as the most distributed species in several regions in China (Chu et al. 2006(Chu et al. , 2010Pan et al. 2011;Wang et al. 2010). This is because MED is pesticide resistant and has a higher competition for field crops (Sun et al. 2013;Horowitz et al. 2020). Recent surveys in China have shown that MED is the dominant cryptic species of B. tabaci in most regions and coexists with MEAM1 in some regions, while some of the 18 native cryptic species identified in China are still dominant in certain regions (Hu et al. 2018(Hu et al. , 2014(Hu et al. , 2015(Hu et al. , 2011Li et al. 2016;Zheng et al. 2016;Liu et al. 2012;Jiu et al. 2017a).
In-depth knowledge of an insect's ecological niche is essential for the formulation of effective pest management strategies. However, B. tabaci is difficult to manage in agricultural landscapes because of its small body size, wide host range, and wide ecological niche (Perring 2001). Moreover, different cryptic species of B. tabaci pose different risks to agriculture due to their biological traits, insecticide resistance, host range, endosymbionts, and ability to transmit viruses (Luan et al. 2013;Liu et al. 2007;Zang et al. 2006;Xu et al. 2011;Shi et al. 2018;Horowitz et al. 2020;Tang et al. 2018;Guo et al. 2019). Previous studies on the distribution and population dynamics of B. tabaci in China have relied mainly on field surveys (Chu et al. 2007;Jiu et al. 2017b;Hu et al. 2014). However, detailed information of the ecological niche complexity of B. tabaci species complex in China is still lacking, and the roles of ecological and anthropogenic factors in the distribution of the cryptic species remain unknown.
Ecological niche modelling (ENM) is a widely accepted approach in studying the potential distribution patterns of species using specific algorithms (Peterson and Soberon 2012). Among ENMs, maximum entropy modelling (Max-Ent) is a general-purpose machine learning method based on the principle of maximum entropy (Phillips et al. 2006). MaxEnt requires only the locations of known occurrences for focus species to express the suitability of each grid cell as a function of the environmental variables at that grid cell, with a high value indicating that the grid cell is predicted to have suitable conditions for the species (Phillips et al. 2006). Due to its easy operation, short model running time, small sample size requirements, and superior performance (Phillips et al. 2006;Hernandez et al. 2006;Pearson et al. 2007;Phillips and Dudik 2008), this method has been widely used to predict the potential distribution of various species, especially invasive organisms, as well as the interactions between species distribution and environments (Orsted and Orsted 2019;Kumar et al. 2014).
Earlier studies have used the MaxEnt model to make effective predictions of the distribution of B. tabaci in North India (Prabhulinga et al. 2017), sub-Saharan Africa (Mudereri et al. 2020), South America (Paredes-Montero et al. 2020), and globally (Ramos et al. 2018). Zhao et al. (2019) used the MaxEnt model combined with GIS technology to predict suitable areas of B. tabaci in China. However, the major drawback of the above studies was that B. tabaci was considered a species rather than a species complex, or only specific cryptic species were examined. For a detailed study of the B. tabaci species complex, it is necessary to consider the different cryptic species.
The ability of an insect pest to survive in an area relies on several factors. Based on the biological characteristics of B. tabaci, the key environmental factors that determine its geographical distribution are temperature and precipitation (Naranjo et al. 2010). Moreover, these variables are correlated with elevation, which may thus influence the occurrence and reproduction of B. tabaci by affecting temperature, precipitation, vegetation, and solar radiation (Azrag et al. 2018;Mudereri et al. 2020). Additionally, since the survival of B. tabaci depends on the host plants, which include several agricultural cash crops, this pest is sensitive to land-use/ land-cover types. Therefore, it is necessary to consider the effects of bioclimatic variables, elevation, and land-use/landcover variables on the distribution and ecological niche of B. tabaci.
In the present study, based on the systematic collection of B. tabaci species complex across China during field surveys and records in the literature, the habitat suitability and distribution of the invasive and native cryptic species of B. tabaci was estimated using MaxEnt software. This study aimed to identify the different ecological niches of B. tabaci cryptic species and the relationship between their presence and key environmental variables.

Sample collection and occurrence data
A large-scale sampling of B. tabaci was conducted in open-fields across China during the peak crop-growing season from 2016 to 2018. We collected whitefly adults from 358 geographical locations in 142 districts in 31 provinces in China. At least 20 individuals were collected from different leaves of separate plants at each site, and at least three sites were surveyed per field. Every individual was identified as a cryptic species using molecular biological methods (Dinsdale et al. 2010;De Barro et al. 2011). We examined published data in the literature on the distribution of B. tabaci in China for sampling information in the wild, including coordinates and cryptic species types, to complement our sampling effort (Hu et al. 2018).
Although extensive investigations and collections have been made, some cryptic species have fewer occurrence records in China. Therefore, four cryptic species with relatively large open-field occurrence records, two invasive cryptic species (MED and MEAM1) and two native cryptic species (China1 and Asia1), were selected for further analysis (Fig. 1). To reduce or minimise sampling deviation and spatial autocorrelation, the spatially rarefy occurrence data tool in the SDM toolbox v2.4 (http:// sdmto olbox. org/) (Brown 2014) was used to rarefy occurrence records in ArcGIS 10.6 (Esri, Redlands, California, USA). After filtering, 234, 59, 38, and 48 presence records of MED, MEAM1, Asia1, and China1, respectively, remained as the occurrence data for habitat suitability modelling (Table S1).

Environmental variables
Three types of environmental variables that were in strong correlation with the occurrence and reproduction of B. tabaci were introduced in this study, including 19 bioclimatic variables (bio1-bio19), elevation (elev), and land-use/land-cover variables (LULC) ( Table 1). The 19 bioclimatic variables and elevation were obtained from the WorldClim Global Climate Database (version 2.1, http:// www. world clim. org) with a spatial resolution of 30 arc-seconds (approximately 1 km at the equator) (Fick and Hijmans 2017). The data of land-use/ land-cover variables were obtained from the Resource and Fig. 1 The occurrence records of the Bemisia. tabaci species complex, including four cryptic species: a MED, b MEAM1, c Asia1, and d China1, which were used as the presence data in MaxEnt modelling 1 3 Environment Science and Data Center (http:// www. resdc. cn), which developed multiple levels for the classification of land-use/land-cover types. The most general or aggregated classification (level I) included broad land-use/land-cover categories, and within each level I class was a number of detailed land-use/land-cover classes (level II) ( Table 2).
When simulating the potential distribution of species, multicollinearity among environmental variables can hinder the analysis of species-environment relationships (Thuiller et al. 2008). To eliminate this problem and improve model performance, the band collection statistics tool in the multivariate toolset from the ArcToolbox of ArcGIS was used to assess the correlations (Pearson correlation coefficient, r) among bioclimatic variables and elevation. If two variables were highly correlated (|r|≥ 0.85), only one variable was included based on realistic biological relevance to B. tabaci (Table S2). Finally, nine environmental variables, namely annual mean temperature (bio1), mean diurnal range (bio2), isothermality (bio3), annual temperature range (bio7), annual precipitation (bio12), precipitation seasonality (bio15), precipitation of the coldest quarter (bio19), elevation, and LULC, were selected for further study (Table 1).

Model development and evaluation
MaxEnt software 3.4.4 (http:// biodi versi tyinf ormat ics. amnh. org/ open_ source/ maxent/) (Phillips et al. 2006) was used to predict the suitable areas and examine the impacts of environmental variables on the distribution of the four cryptic species of B. tabaci in China. This program has the advantages mentioned above, as well as being most suited to our research which had only presence data available for B. tabaci. The models were set to 10 replicates with subsample method, and each run implemented 75% of the occurrence points as the training data and the remaining 25% as the test data. Response curves were created to show how the relative probability of occurrence depended on the value of each environmental variable. The jackknife technique was used to measure the relative importance of a particular variable by comparing the gain of the model using a single variable, excluding the variable and all the variables (Phillips et al. 2006). The maximum number of iterations was increased from 500 (default) to 5000 to optimise convergence. The cloglog output format was checked for model analysis.
In the MaxEnt algorithm, feature classes (FCs) and regularisation multiplier (RM) are two important settings to adjust as, depending on the dataset, their values can greatly affect model performance (Muscarella et al. 2014). To achieve models with the greatest predictive power, the R package "ENMeval" (Muscarella et al. 2014) was used to select the optimal parameter combinations of FCs and RM for each cryptic species. MaxEnt provides five feature types: linear (L), quadratic (Q), hinge (H), product (P), and threshold (T) (Phillips and Dudik 2008). Nine different FCs Mean diurnal range (mean of monthly (max temp-min temp)) bio3 Isothermality (bio2/bio7) (× 100) bio4 Temperature seasonality (standard deviation × 100) bio5 Maximum temperature of the warmest month bio6 Minimum temperature of the coldest month bio7 Annual temperature range (bio5-bio6) bio8 Mean temperature of the wettest quarter bio9 Mean temperature of the driest quarter bio10 Mean temperature of the warmest quarter bio11 Mean temperature of the coldest quarter bio12 Annual precipitation bio13 Precipitation of the wettest month bio14 Precipitation of the driest month bio15 Precipitation seasonality (coefficient of variation) bio16 Precipitation of the wettest quarter bio17 Precipitation of the driest quarter bio18 Precipitation of the warmest quarter bio19 Precipitation of the coldest quarter Elevation elev Ground height above sea level Land-use/land-cover variables LULC Land-use/land-cover types (see Table 2 for additional details) (L, LQ, LH, LQH, LQP, QHP, LQHP, QHPT, and LQPHT) were tested with RM ranging from 0.5 to 4.0 in increments of 0.5. The "Checkerboard2" approach was applied by calculating the lowest delta value of the Akaike information criterion coefficient (AICc) to select the final MaxEnt models for each cryptic species. Finally, four sets of optimal parameter combinations of FCs and RM, namely LQP with 1.5, L with 0.5, LQP with 4.0, and LQHP with 2.0, were selected to execute the final model for MED, MEAM1, Asia1, and China1, respectively (Fig. 2). The MaxEnt models were evaluated using a threshold independent receiver operating characteristic (ROC) curve, and the area under the curve (AUC) was positively correlated with the performance of the prediction model (Peterson et al. 2008). Thus, the model performance was classified depending on the AUC value as failing (0.5-0.6), poor (0.6-0.7), fair (0.7-0.8), good (0.8-0.9), or excellent (0.9-1.0) (Phillips et al. 2006).

Habitat and ecological niche analysis
The maximum test sensitivity plus specificity (MTSPS) threshold was chosen to determine the suitable and unsuitable habitats of B. tabaci, as the thresholds maximising sensitivity and specificity performed well on presence-only datasets ). Based on this, four suitability levels of distribution regions were defined for each cryptic species: (1) unsuitable (P < MTSPS), (2) low suitable (MTSPS ≤ P < 0.5), (3) moderately suitable (0.5 ≤ P < 0.7), and (4) highly suitable (P ≥ 0.7) (Ramos et al. 2019).
To compare the differences in the geographical distribution of B. tabaci cryptic species, the MaxEnt outputs were imported into ArcGIS for suitability classification and visual interpretation, and the suitable habitat areas were further counted using the zonal statistics tool. Subsequently, their ecological niche differences, including niche breadth, niche overlap, and range overlap tests, were performed using ENMTools 1.3 (Warren et al. 2010) based on MaxEnt outputs. In the niche overlap test, Schoener's D index (Schoener 1968) and Hellinger's-based I index (Warren et al. 2008) were used to measure the similarities between each pair of cryptic species habitat suitability models. These values were between 0 and 1, with higher values corresponding to higher indicators (Warren et al. 2010).

Quantifying importance of environmental variables
The dominant environmental variables affecting the distribution of the four cryptic species were analysed based on relevant measures (i.e. jackknife results, percentage contribution, and permutation importance), and the relationships between the probability of their presence and these variables were described from the response curves. As the parameters (FCs and RM) may significantly impact the response curves, four supplementary experiments were established. In different experiment, one of the optimal parameter combinations mentioned in Sect. 2.3 was applied as fixed settings to derive the response curves of environmental variables for the four cryptic species. In addition, the adaptation ranges of the nine environmental variables for the four cryptic species were extracted and then calculated separately using the extract by points tool in ArcGIS.

Model performance
The MaxEnt models had high accuracy with average AUC values higher than 0.9, MED (0.922), MEAM1 (0.924), Asia1 (0.971), and China1 (0.953) (Fig. 3). The habitat distributions of the four cryptic species generated using Max-Ent modelling were reliable and could be used for further analysis.

Suitable habitat distribution
The models showed areas of low to high habitat suitability for the four cryptic species of B. tabaci in China (Fig. 4). Generally, the distribution area of the invasive cryptic species was larger than that of the native cryptic species. Among the four cryptic species, the two invasive cryptic species had similar habitat distribution. Specifically, MED and MEAM1 were distributed almost nationwide in China, except for parts of northwest, southwest, and northeast China, with the moderately and highly suitable habitats mainly located in northern, central, and eastern China for MED, and in northern, central, eastern, and southern China for MEAM1 ( Fig. 4a-b). Asia1 was mainly distributed in southwest and southern China (Fig. 4c), whereas China1 was mainly distributed in eastern, central, and southern China (Fig. 4d). In terms of the areas of suitable habitat, MED was the highest among the four cryptic species, followed by MEAM1, China1, and Asia1 in descending order (Table 3).

Comparisons of ecological niche
Ecological niche analysis showed that MED had the highest niche breadth followed by the remaining three in the order of MEAM1 > China1 > Asia1 (Table 4).
We observed some degree of niche overlap between the four cryptic species. Both the D and I indices indicate a high degree of niche overlap between MED and MEAM1 Fig. 2 The results of ENMeval package for the four cryptic species of Bemisia. tabaci, which showed that the optimal parameter combinations of FCs and RM in the MaxEnt models were: a LQP and 1.5 for MED, b L and 0.5 for MEAM1, c LQP and 4.0 for Asia1, and d LQHP and 2.0 for China1, respectively (D > 0.7, I > 0.9), indicating that their ecological niches might be similar, but not identical, followed by MEAM1 and China1 (D > 0.5, I > 0.8). However, the other pairs of cryptic species did not share similar niches (Table 4).
Range coverage analysis showed that MED and MEAM1 had the highest range coverage, followed by MED and Asia1, MEAM1 and Asia1, MEAM1 and China1, MED and China1, and Asia1 and China1, in descending order (Table 4). There was a large niche overlap between Asia1 and China1, while the range overlap was small (Table 4).

Important environmental variables
Based on the results of the modelling, the important environmental variables affecting the distribution of the cryptic species were different. However, the important environment variables affecting the distribution of MED and MEAM1 were similar. LULC, elevation, and bio1 were the three most important factors among the nine environmental variables, with a cumulative contribution of 70.7% for MED and 89.4% for MEAM1 ( Fig. 5a-b, Table 5). The permutation analysis ranked these variables as the three most important variables, although in a different order for MED and MEAM1 (Table 5). Furthermore, the result of the jackknife test shows that LULC contained the most useful information of a single variable and may contribute to additional information that is not present in the other variables ( Fig. 5a-b). However, the results for the two native cryptic species were different from those obtained for the invasive cryptic species. Bio1 and bio3 were the main factors affecting the distribution of Asia1, with a cumulative contribution of 94.0% (Table 5). Additionally, bio1 had the highest permutation importance for Asia1 (Table 5) and contained the most useful information compared with the other variables (Fig. 5c). For China1, the main factors were bio19, elevation, and bio12, with a cumulative contribution of 73.5% (Fig. 5d, Table 5). Among the variables, bio1 and elevation had a higher permutation importance ( Table 5). The result of the jackknife test shows that bio1 had the most useful information as a single variable, and bio7 had the most information that was not present in the other variables (Fig. 5d). Overall, temperature-related variables contributed more than precipitation-related variables to the models of the four cryptic species.

Response to environmental variables
The trends of response curves show the relationships between the probability of presence of the cryptic species     . 6 Response curves of the environmental variables most related to the distribution of MED in elevation, dropping to almost 0 at approximately above 3000 m (Fig. 6b), but exhibited a quadratic pattern with an increase in bio1, peaking at approximately 21.5 ℃ and declining at higher temperatures (Fig. 6c). The probability of the presence of MEAM1 shows the same response to LULC (Fig. 7a) and elevation (Fig. 7c) as MED and increased with an increase in bio1 (Fig. 7b). The probability of the presence of Asia1 increased with an increase in bio1 (Fig. 8a) and bio3 (Fig. 8b), and the most important land types, with the exception of water, were the same with MED and MEAM1, although there were differences in the response curves of LULC (Fig. 8c). The probability of the presence of China1 increased with an increase in bio19 (Fig. 9a), but decreased with an increase in elevation, dropping to almost 0 at above approximately 2000 m (Fig. 9b). In terms of bio12, China1 showed a quadratic pattern with an increase in this variable, peaking at approximately 1300 mm, and declining with higher precipitation (Fig. 9c).
In four groups of supplementary models with fixed parameter combinations, the different trends of the response curves of the same variable indicate the different environmental requirements of the different cryptic species (Fig.  S1-S4). The results show that the parameter combination (FCs and RM) did affect the response curves; however, the response of the same variable for a particular cryptic species remained the same in most groups, such as bio1, bio2, and bio7. Specifically, the responses to elevation for the four cryptic species remained almost identical. The trends of  response curves of bio3 for the two native cryptic species were always opposite, while that of bio12 for MEAM1 and Asia1, and bio15 for MEAM1 and China were, respectively, consistent.

Adaptation for environmental variables
The adaptation ranges of the main environmental variables for the four cryptic species further illustrate the niche complexity of the B. tabaci species complex. There were significant differences in the adaptation range of the nine environmental variables (Fig. 10). The difference in the range of adaptation of the four cryptic species to bio1 indicates that invasive cryptic species were better adapted to lower temperatures than native cryptic species (Fig. 10a). The difference in bio2, bio3, and bio7 indicates that Asia1 had the lowest adaption to temperature variations (Fig. 10b-d). For precipitation-related variables, the difference in adaptation range for bio12 and bio15 revealed that MED had a higher adaptation to areas with low precipitation (Fig. 10e) and large seasonal variations in precipitation than the other cryptic species (Fig. 10f). China1 did not seem to be able to adapt to drought in the coldest quarter, judging from the result of bio19 (Fig. 10g). The result for elevation indicates that Asia1 had a wider adaptation range to elevation than the other three cryptic species, while China1 was only suitable for survival at lower elevation (Fig. 10h). MED had the highest adaptation to a wide range of LULC, which were mainly agriculture (Level I number 1) and urban or builtup (Level I number 5). However, the other cryptic species had similar land preferences, which were mainly urban or built-up (Level I number 5) for MEAM1, and agriculture (Level I number 1) for Asia1 and China1 (Fig. 10i).

Discussion
In this study, we sought to create ecological niche models for different cryptic species of B. tabaci in China. Their occurrence data were collected mainly through extensive field sampling across China and further identification of cryptic Fig. 10 The adaptation range of different environmental variables for the four cryptic species of Bemisia. tabaci in China, a annual mean temperature (bio1), b mean diurnal range (bio2), c isothermality (bio3), d annual temperature range (bio7), e annual precipitation (bio12), f precipitation seasonality (bio15), g precipitation of coldest quarter (bio19), h elevation (elev), and i land-use/land-cover types (LULC) species using reliable molecular means. However, the occurrence data are often clustered in geographical space, and models created using such data may over-fit towards environmental biases and model performance values may be inflated, producing erroneous results (Veloz 2009;Boria et al. 2014). A professional tool was used to address this issue by spatially filtering the locality data to reduce occurrence localities to a single point within the specified distance, ensuring that only one point per cell was available (Brown 2014). Moreover, apart from bioclimatic variables, elevation and land-use/land-cover variables with high spatial resolution were added to the MaxEnt model to improve the accuracy of the prediction of distribution for the cryptic species. Hence, combining the AUC values, the models used in the present study had a high degree of reliability.
Among the cryptic species of B. tabaci, MED and MEAM1 are the most invasive and destructive, affecting agricultural production globally (De . During field surveys, we found that the two invasive cryptic species are distributed in most parts of China and the distribution of MED has gradually exceeded that of MEAM1 in recent years (Jiu et al. 2017b;Hu et al. 2018;Tang et al. 2020). However, some native cryptic species are still dominant in some areas (Hu et al. 2015;Tang et al. 2018). The modelling results suggest that the four cryptic species examined in the present study, two invasive cryptic species (MED and MEAM1) and two native cryptic species (China1 and Asia1), had specific habitat distribution patterns and ecological niches in China, which are consistent with the results of field surveys. Our study demonstrate that there are differences in the distribution patterns and ecological niches of invasive and native cryptic species of B. tabaci.
This study provides the first predictive distribution maps for different cryptic species of B. tabaci. The two invasive cryptic species have the widest distributions patterns, covering almost all of China except for the Qinghai-Tibetan Plateau and northeast extremely cold area, while their high suitability areas are not overlapping. In contrast, the native cryptic species have smaller distribution ranges, each concentrated in certain areas of China. Moreover, based on the field observations and model results of this study, the coexistence of B. tabaci cryptic species is a common condition, although invasions may result in the replacement of native species by alien-species (Guo et al. 2021). Niche partitioning and stochastic processes across landscapes can promote this process (Crowder et al. 2011).
The niche breadth of the invasive cryptic species significantly exceeded that of the native cryptic species, suggesting that the invasive species were better adapted to different environments, which may be a reason why the former is more threatening to crop production than the latter. The niche breadth of MED was larger than that of MEAM1, which may be a reason why MED has gradually replaced MEAM1 in some regions of China in recent years. Furthermore, the two invasive cryptic species have greater niche overlap and range overlap, indicating that there may be a relatively intense competitive ecological niche relationship between them. However, the environmental and ecological factors behind MED replacement of MEAM1, as found in some of the field surveys, need to be further investigated.
The four cryptic species had specific sets of dominant environmental variables and exhibited different adaptations to these variables, which may explain their differences in distribution and ecological niche. The type of land-use/landcover is an important factor affecting the occurrence of MED and MEAM1 in China, with urban or built-up land having the highest suitability for the survival of the pest, followed by agricultural land. This indicates that the survival of the invasive cryptic species of B. tabaci depends on the host plant, which may be closely related to human activities. The results of the present study show that B. tabaci occur mainly at low elevations and its probability of survival decreases with an increase in elevation, which is consistent with a previous study in China (Zhao et al. 2019). However, some cryptic species may have a wider adaptability to changes in elevation. Among the bioclimatic variables, temperaturerelated variables had a greater effect on the distribution of B. tabaci than precipitation-related variables, with precipitation barely playing a role in the distribution of Asia1. Studies have shown that MED and MEAM1 can survive in the tropics where the temperature is higher (Ramos et al. 2018;Xiao et al. 2016). The results of the present study indicated that low temperature is a limiting factor for the occurrence and survival of MED and MEAM1 in the wild in China. However, with climate change causing an increase in global and regional temperatures, MED and MEAM1 may be able to survive in areas previously not suitable for them. Thus, how the distribution ranges of B. tabaci cryptic species will change under future climate change conditions deserves further study.
Our results provide insight to the ecological characteristics of the B. tabaci species complex, suggesting that future studies on this pest need to consider the cryptic species in isolation, especially in terms of ecology. Moreover, prevention and field management should be strengthened to combat the spread of invasive cryptic species into suitable areas where they do not yet occur, particularly at sites under high risk of the invasive cryptic species. Meanwhile, our results also provide insight for different regional management strategies in China as different cryptic species showed different resistance to some insecticides and sensitivity to predominant parasitoids (Dennehy et al. 2010;Zhang et al. 2015;Horowitz et al. 2020).

Conclusion
The present study provides estimates of the potential suitability of invasive and native cryptic species of B. tabaci in China using MaxEnt and further identified the differences in their ecological niches. Our analysis revealed that the B. tabaci species complex exhibits significant differences in suitable area distribution, important environmental variables affecting their distribution, adaptability to environmental variables, ecological niche breadth, and niche overlap. Specifically, the invasive cryptic species are distributed over most of China, while the native cryptic species are limited to certain areas of China. The distribution range and niche breadth of the invasive cryptic species excess that of the native cryptic species in the order of MED > MEAM1 > China1 > Asia1, and there are different degrees of niche overlap and range overlap among them. These results might explain the phenomenon found in field surveys from the perspective of the ecological niches and support that the B. tabaci species complex occupies a complex ecological niche in China. The findings of this study improve our understanding of the ecological characteristics of B. tabaci species complex, which will be useful in the development of prevention and control strategies for this pest in China. In addition, this study also provides an approach for conducting related research in other regions of the world.

Author contributions
Y.X., C.L., and L.J. conceived and designed research; Y.X., Y.W, and Y.Z. collected and tested the samples; Y.X. conducted experiments; Y.X. and C.L. analysed data; Y.X. wrote the manuscript. All authors read and approved the manuscript.