Adaptation of acaricide stress facilitates Tetranychus urticae expanding against Tetranychus cinnabarinus in China

Abstract The two‐spotted spider mite, Tetranychus urticae, and the carmine spider mite, Tetranychus cinnabarinus, are invasive and native species in China, respectively. Compared with T. cinnabarinus, T. urticae has expanded into most parts of China and has become the dominant species of spider mite since 1983, when it was first reported in China. However, the mechanism of the demographic conversion has not been illuminated. In this study, one T. urticae field population and one T. cinnabarinus field population were isolated from the same plant in the same field, and the toxicological characteristics were compared between these two species. Laboratory bioassays demonstrated that T. urticae was more tolerant to commonly used acaricides than T. cinnabarinus. The activities of detoxification enzymes were significantly greater in T. urticae, and the fold changes of enzymes activities in T. urticae were also greater following exposure to acaricides. Furthermore, more metabolism‐related genes were upregulated at a basal level, and more genes were induced in T. urticae following exposure to acaricides. The comparison of proteins and genes between both species led credence to the hypothesis that T. urticae was more resistant to acaricides, which was the reason explaining the expansion of invasive T. urticae against native T. cinnabarinus. Laboratory simulation experiments demonstrated that following the application of acaricides, the composition of a mixed T. urticae/T. cinnabarinus population would change from a T. cinnabarinus‐dominant to a T. urticae‐dominant population. This study not only reveals that T. urticae possesses stronger detoxification capacity than its sibling species T. cinnabarinus, which facilitated its persistent expansion in China, but also points to the need to accurately identify Tetranychus species and to develop species‐specific management strategies for these pests.

In China, the two forms coexist, where the red form T. cinnabarinus is distributed throughout China and is considered to be native, while the green form T. urticae, which was first reported in 1983 in Beijing, is considered to be invasive (Dong, Guo, & Niu, 1987;Sun, Lian, Navajas, & Hong, 2012). Tetranychus urticae has recently expanded its distribution from its putative area of introduction in Beijing to many parts of the country, including Hebei, Liaoning, Jilin, Gansu, Anhui, and Yunnan provinces and elsewhere (Meng, Wang, Jiang, & Yi, 2001;Sun et al., 2012).
Unexpectedly, the invading T. urticae has become the most important mite pest in apple orchards in the north of China (Cai, Cheng, & Sha, 2003). A recent study by Wang indicated that T. urticae and T. truncatus were the dominant species on vegetables in some areas in Beijing and Hebei, revising the traditional opinion that T. cinnabarinus was the major species on common vegetable plants (Wang, Zhang, Wu, Xie, & Xu, 2013).
Most toxicity studies have concluded that there are differential responses of the two mites to commonly used acaricides (including plant resource pesticides). For example, T. urticae has greater resistance to abamectin, spiromesifen, etoxazole, hexythiazox on strawberries than T. cinnabarinus in California, USA (Bi, Niu, Yu, & Toscano, 2016).
T. urticae was also identified to be more tolerant than T. cinnabarinus to traditional acaricides, such as fenazaquin, pyridaben, propargite, azocyclotin, and hexythiazox (Gu, Zhang, Zhao, & Li, 2000). In comparison with T. cinnabarinus, T. urticae has a significantly decreased susceptibility to six acaricides, which has been shown to be true in both slidedip immersion and leaf-dipping method testing, where remarkably, the tolerance of T. urticae to abamectin was 2,575-fold higher than that of T. cinnabarinus (Zhao, Zhou, & Ren, 2006). Tetranychus urticae has also been shown to be more tolerant to natural products, where the LC 50 values for crude extracts from Wikstroemia chamaedaphne using chloroform, petroleum ether, or ethyl acetate were higher for T. urticae than for T. cinnabarinus (You-Nian et al., 2010).
Metabolic resistance has been reported worldwide and usually involves detoxification enzymes such as cytochrome P450 monooxygenases (P450s or CYPs for genes), carboxy/cholinesterases (CarEs or CCEs), and glutathione S-transferases (GSTs) (Van Leeuwen et al., 2010;Xu et al., 2014). Gene expression analysis is widely used to reveal regulatory mechanisms that control cellular processes in animal, plants, and microbes (Van Leeuwen, Dermauw, Grbic, Tirry, & Feyereisen, 2013), including the elucidation of the gene expression profiles of detoxification genes involved in the metabolic resistance of insecticides (Strode et al., 2008). In particular, recently developed massively parallel RNA-Seq deep sequencing and digital gene expression (DGE) testing have substantially changed the way resistance-relevant genes in insects are identified and characterized because these methods facilitate the investigation of the functional complexity of the transcriptome (Metzker, 2010). Moreover, RNA-Seq offers a great depth of sequence coverage with reduced variability (Zhou et al., 2010).
Following the initial discovery of T. urticae, the range of this species has expanded throughout China and has become the dominant spider mite. The causes of the demographic change are not always clear and may be attributed to various biological mechanisms or anthropogenic factors, including differences in their respective reproductive success, and susceptibility to acaricides. Similar effects can be seen in other insect species, where it has been reported that the application of insecticides can lead to rapid shifts in the composition of leafminer complexes in laboratory and field-based experiments (Gao, Reitz, Wei, Yu, & Lei, 2012). Although the greater resistance of T. urticae versus T. cinnabarinus to acaricides is generally recognized, it is unclear whether acaricides facilitate the expansion of T. urticae, resulting in T. urticae being the dominant species of spider mites.
In our study, we compare the relative fitness at different temperatures between T. urticae and T. cinnabarinus, as well as the susceptibility of the two mites collected from the same field crop to commonly used acaricides. Furthermore, the influence of acaricides on changes in the population composition of mixed populations of T. urticae and T. cinnabarinus was simulated under laboratory conditions. Finally, the biochemical characteristics and gene expression profiles of metabolic detoxification enzymes in both mite species were characterized and compared. The goal of our study is to elucidate the mechanism(s) that has resulted in the continuous expansion of T. urticae throughout China to become the dominant pest mite species.

| Mites and acaricides
Tetranychus urticae population and Tetranychus cinnabarinus population were originally collected from the same commercial rose field in Kunming, Yunnan, China (25°23′N, 102°42′E) in May 2014, and named Tu-YN, Tc-YN, respectively ( Figure 1). Yunnan is located within a subtropical zone with a yearly average temperature of 23°C and maximum temperature of 38.5°C. The two spider mite populations were regarded as having the same acaricide-exposure background as they were collected from the same crop of roses, which had been sprayed mainly with abamectin and occasionally with other acaricides (e.g., propargite, pyridaben, or cyflumetofen). The Tc-SS strain was used as a reference colony of T. cinnabarinus, which was originally established from more than 1,500 mites collected from a cowpea field in Beibei, Chongqing, China, in 1998 (Li, Chen et al., 2009;Li, Lu et al., 2009), and has since been continuously reared in the laboratory. All strains and stock cultures were reared on detached cowpea (Vigna unguiculata L.) leaves placed on water-soaked cotton in petri dishes (9 cm diameter) in growth chambers.

| Number of setae on tibia of leg I in female
Mites were slide-mounted with Hoyer's solution and cured on a Flatting table (Leica M205A, Germany) at 50°C for 4-5 days. Mounted mites were observed with an Olympus BX51 (Japan) equipped with differential interference contrast. The tibia setation follows Kuang and Cheng (1990) with modifications of Norton (Norton, 1998) concerning coxal setae. Thirty adult females picked randomly were observed from the corresponding populations (Tc-YN, Tu-YN, and Tc-SS, respectively), which had been reared under the same condition for 12 months.

| Hybridization
Tetranychus spp. have haplodiploid sex determination, where unfertilized eggs develop into males, while zygotes will develop into females. Therefore, reproductive isolation can be detected based on the sex ratio of the offspring, as females can only be produced through a viable cross. Cross-mating of T. urticae to T. cinnabarinus was conducted in six-cm petri dish arenas containing a bean leaf disk (4 cm diameter) on water-soaked cotton. One female deutonymph (Tc-YN, Tu-YN, or Tc-SS) was transferred to each arena, and then three virgin males from another mite species were added. The males and female were removed 2 and 5 d after adult eclosion, respectively. The numbers of F 1 females and males was recorded and if adult females were present among the F 1 offspring, they were allowed to in-cross and the numbers of F 2 females and males were recorded through to the F 3 generation.

| Data collection of the expansion of Tetranychus urticae in China
Studies regarding T. urticae and T. cinnabarinus in China were collected by retrieving the term "Tetranychus urticae" and/or "Tetranychus cinnabarinus" in the CNKI database in Chinese and in the SCI database.
Based on the historical reports from the literature, the geographical distributions of T. cinnabarinus and T. urticae were determined based on the reported locations and dates from which the mites were collected. The map of distribution was created using ArcGIS9.3 software (Environmental Systems Research institute, Redlands, CA, USA) (http://www.arcgis.com) and Adobe Illustrator CS5 (v 15.1.0.39) software (http://www.adobe.com/products/illustrator.html).

| Biology of Tc-YN and Tu-YN under the different temperature conditions
Life-history variables of Tc-YN and Tu-YN were measured using arenas, each consisting of a bean leaf disk (4 cm diameter) on watersoaked cotton in a petri dish (6 cm diameter). The arenas were kept in climatic chambers at 16 ± 1°C/26 ± 1°C/33 ± 1°C, respectively, 50 ± 5% RH, and a 14:10 (L:D) regime. Five adult females were transferred to each arena for 4 hr to lay eggs, after which the females were removed and a single egg was left in each arena. A total of 50 arenas were prepared for each species per temperature. The arenas were examined every 24 hr, and the duration and mortality for each development stage were recorded.
F I G U R E 1 The two-spotted spider mite, Tetranychus urticae, and the carmine spider mite, Tetranychus cinnabarinus are both economically significant species of the genus Tetranychus, belonging to the class Arachnida, infraclass Acari, order Prostigmata, and family Tetranychidae Life-tables were prepared from the estimated rates for age-specific survival (l x ) and age-specific fecundity (m x ), where l x and m x are the proportion of surviving females at age x and the number of oviposition events produced per female in the age interval x, respectively. The sum of the products (l x m x ) then provides the net reproductive rate (R 0 ). The mean generation time (T, days) was calculated as: T = Σl x m x x/R 0 .
The innate capacity for increase (r m ) was calculated as: r m = lnR 0 /T (where e is the natural logarithm base), and the population doubling time (t, days) was computed using the formula: t = ln2/r m . Relative fitness (R f ) was estimated as follows: R f = (R 0 of Tu-YN)/(R 0 of Tc-YN), an R f > 1 suggests that the Tu-YN has a fitness advantage, whereas an R f < 1 suggests that Tu-YN has a fitness defect compared with Tc-YN (Groters, Tabashnik, & Finson, 1994;Li, Gao, Zheng, & Liang, 2000).

| Toxicological tests and acaricide stimulation
Median lethal concentration (LC 50 ) values for adult mites were measured using the modified residual coated vial (RCV) method recommended by Van Leeuwen, Stillatus, and Tirry (2004). Five acaricides (abamectin, tebufenpyrad, cyflumetofen, propargite, and bifenazate), which could kill the larva, were selected to carry out the larval bioassays, and these tests were based on the method described by Knight, Beers, Hoyt, and Riedl (1990). Thirty larvae were placed on a bean leaf disk (4 cm diameter) on water-soaked cotton in a petri dish (6 cm diameter). Prepared solutions of acaricides (1 ml) were sprayed onto the leaf disks containing the mites using a Potter Spray Tower (Burkard Manufacturing, Rickmansworth, Herts, UK) at 68.9 kPa. The acaricide solutions were diluted in distilled water to at least seven concentrations that resulted in mite mortality ranging from 20 to 85%. Distilled water alone was used as a control, and all doses were independently replicated three times. Following acaricide application, mites were maintained on the leaf disks at 26 ± 1°C and a 14:10 (L:D) regime. Individual mites were assessed for mortality under a dissecting microscope 24-48 hr after treatment, when all larvae in the control treatment had developed into nymphs. LC 50 values and their 95% confidence limits were calculated from probit regression using probit analysis software (POLO, leOra Software, Berkeley, CA, USA). If the 95% confidence limits of the LC 50 values did not overlap between Tc-YN and Tu-YN, then the differential acaricide susceptibility was considered to be significant.
To test the response of mites to abamectin-, fenpropathrin-, and tebufenpyrad stimulation, newly emerged adults (2-3 days old) Tu-YN and Tc-YN were treated for 6 hr with an LC 30 dosage as described by Shen, Shi, Xu, and He (2014). For this, adults were exposed to tubes coated with the acaricide and acetone served as the control. Surviving mites were collected for enzyme sources and RNA extraction. Each sample was independently replicated three times for enzyme activity assay and RNA sequencing.

| Mixed population experiments
To generate mixed species populations, five male and five female adults (newly eclosed virgins) from each of the Tu-YN and Tc-YN populations were added to the same bean leaf disk (6 cm diameter) on water-soaked cotton in a petri dish (9 cm diameter). Mites were held on the leaf disk at 26 ± 1°C and a 14:10 (L:D) photoperiod and allowed to mate and expand their populations for 15 days. After 15 days, the numbers of female adults belong to each species were recorded and the leaf disks were sprayed with either abamectin (4.5 mg/L, the recommended field dosage), cyflumetofen (80 mg/L, the recommended field dosage), or distilled water (as a control) using a Potter Spray Tower (Burkard Manufacturing, Rickmansworth, Herts, UK) at 68.9 kPa. The spray applications were repeated in intervals of 15 days for a total three times. The numbers of female adult mites of each species were measured immediately prior to each spray, and again for a final reading at 15 days following the last spray. All treatments were independently replicated three times.
The supernatant was then used for testing using nitroanisole (0.05 M in acetone) as the substrate, and NADPH was added to reaction. The reaction was allowed to proceed for 30 min at 37°C and then stopped with 1 M hydrochloric acid, extracted over chloroform and neutralized with 0.5 M NaOH. The optical densities of the reactions were then measured at 400 nm using a microplate reader (TECAN Co.). Wells without extracted enzyme were used as controls, and the amount of protein in the enzyme source was determined using the Bradford method, with bovine serum albumin as the standard. The specific activity was determined according to a nitrophenol standard curve and the protein concentration of the sample.

| Glutathione S-transferase activity assay
Glutathione S-transferases activity was determined according to the method of Habig, Pabst, and Jakoby (1974) with modifications by Stumpf & Nauen (2002). A total of 200 female adult mites (Tu-YN, Tc-YN) were homogenized on ice in 1.5 ml PBS (0.04 M, pH 6.5), followed by centrifugation at 10,000 × g for 10 min at 4°C. CDNB (0.6 mM) and GSH (6 M) were used as substrates and added to the supernatant, which was incubated for 20 min at 37°C during which time, GSTs reacted with the reduced GSH. The optical density at 340 nm was immediately recorded at intervals of 30 s for 5 min using a microplate reader. The results were determined based on the protein concentration of sample, and the specific activity was converted from the OD value.

| Carboxylesterase activity assay
The method reported by Van Asperen (Van Asperen, 1962) was adopted for testing CarE activity, in which 200 female adult mites (Tu-YN, Tc-YN) were homogenized on ice in 1.5 ml PBS (0.04 M, pH 7.0), then centrifuged at 10,000 × g for 10 min at 4°C. The supernatant was removed and placed on ice, and α-naphthyl acetate (3 × 10 −4 M) and 10 −4 M physostigmine were used as the substrates. The reaction was incubated for 10 min at 37°C, then the color-developing agent was added (mixed as follows: mass fraction 5% SDS: mass fraction 1% fast blue B salt = 5:2 (v/v)), and the OD value at 600 nm was recorded immediately. The specific activity of CarE was calculated based on the α-naphthol standard curve and the protein concentration of sample.
Genomic DNA was removed using a genomic DNA elimination column supplied with the kit. The quality of the RNA sample was verified by ensuring that the OD 260/280 was within the range of 1.8-2.2 as measured by the NanoVue UV-vis spectrophotometer (GE Healthcare Bio-Science, Uppsala, Sweden), and qualified samples also had a 28S to 18S rRNA ratio above 1.0 as measured by a 1% agarose gel electrophoresis and the Agilent 2100 Bioanalyzer (Palo Alto, CA, USA) with a minimum integrity value of 8.

| cDNA library preparation and RNA-Seq
Total RNA pools from each sample were then used for preparing cDNA libraries using the Ion Total RNA-Seq kit V2 (Life Technologies Corporation, CA, USA). Double-stranded cDNA was ligated to barcoded adapters and was sequenced using the Ion PI ™ Chip (Ion torrent, Life technologies, CA, USA) at the Beijing Genomics Institute (Shenzhen, China). Libraries were run at a concentration of 4-5 pM. To ensure the accuracy of subsequent analysis, raw sequences were cleaned to remove adaptors and sequencing errors. Reads were removed if they contained the sequencing adaptor, more than 5% unknown nucleotides, or more than 20% of bases of low quality. This output was called clean reads, which was used for subsequent downstream analyses.

| Gene mapping
The entire assembled genes were used to search for the best-hit homologous proteins (BLASTX cutoff e-value 1.0E-5) in the T. urticae genome (Grbić et al., 2011). Ortholog prediction was performed by performing BLASTX and TBLASTN bidirectional comparisons between Tc-YN and Tu-YN (with a threshold cutoff e-value ≤1. 0E-5) to identify the hits within the two species. between two different DGE libraries was conducted using the method described by Audic and Claverie (Audic & Claverie, 1997).

| Differential gene expression analysis
The false discovery rate (FDR) method determines the p-value threshold for multiple testing by controlling the FDR value. The criteria of FDR < 0.001 and the absolute value of log 2 ratio >1 were used to judge the significance of gene expression differences. The full dataset has been submitted to NCBI-GEO database under the experiment ID GSE75529. The DEG functional features were analyzed according to cluster of orthologous groups of proteins (COG) and KEGG orthology (E-value, 10E-5). The Blast2GO program (Conesa et al., 2005) was employed to obtain Gene Ontology (GO) annotations for the DEGs.

| Cloning and sequence comparison of target of acaricide genes for Tu-YN and Tc-YN
The coding sequence (CDS) information of 14 genes (VGSC, GluCl-01-05, Rdl1-3, SdhA-E) (Table S10) targeted by acaricides was obtained from the NCBI database for Tc-YN, and the T. urticae genome (Grbić et al., 2011) for Tu-YN and used to design PCR primers (Table S10)

| Statistical analyses
The statistically significant differences of development duration, enzyme activity, and individual numbers in mixed population were calculated using independent-sample t-tests for all two-sample comparisons with a p-value <.05 using SPSS19.0 software, CSPSS, Chicago, IL, USA.

| Species identification
Following field collection, adult female mites with two spots and red bodies were identified as T. cinnabarinus and designated as Tc-YN, while adult female mites with two spots and green bodies were identified as T. urticae and designated as Tu-YN. From these two mite populations, 30 individuals were randomly selected from each for the counting of the setae on the tibia of leg I. In the Tc-YN and Tc-SS, the tibiae of the first pair legs of most adult females had 10 setae, with some individuals having 12 or 13 setae (Table 1). In the Tu-YN population, however, only 10 setae were observed (Table 1, Figure   S1), which is a feature of T. urticae population (Boudreaux, 1956;Kuang & Cheng, 1990). The morphological characters observed were similar to those reported by Kuang and Cheng (1990) and Zhang and Jacobson (2000) who reported that the number of setae on tibia I in adult females was a very useful and convenient method for separating T. cinnabarinus and T. urticae populations.
The results of cross-breeding between T. urticae and T. cinnabarinus were presented in Table 2. Reciprocal cross-mating between two T. cinnabarinus populations originating from the field (Tc-YN)and laboratory (Tc-SS)-produced normal progeny females (F 1 to F 3 ) showing no reproductive isolation between these two populations.
Cross-mating between the two T. cinnabarinus populations and Tu-YN produced only males in the F 1 generation, with the exception of one Tc-SS × Tu-YN' mating, which generated F 1 females with an unusual wax yellow color that could not generate F 2 female offspring (Table 2, Figure S2). The results suggested that reproductive isolation between T. urticae and T. cinnabarinus was complete.

| Comparison of biological characteristics between Tc-YN and Tu-YN
Comparisons between Tc-YN and Tu-YN at 16°C indicated significant differences in the development duration (the duration of the period from egg to adult female), the net reproductive rate (R 0 ), the capacity for increase (r m ), the population doubling time (t) and the average generation time (T). The main contrast between Tc-YN and Tu-YN was observed with respect to R 0 , and the value for this for Tc-YN was estimated to be 30.7% lower than for Tu-YN. The relative fitness (R f ) observed for Tu-YN was 1.44, indicating the presence of a fitness advantage in Tu-YN at lower temperatures compared with Tc-YN (Table 3). At 26°C, no significant differences between Tc-YN and Tu-YN were recorded for all population parameters (  (Table 3).

| Susceptibility of Tc-YN and Tu-YN to acaricides
The toxicities of abamectin and seven other acaricides to Tu-YN and In T. cinnabarinus, the number of the female mite is greater than that of male, and the ratio range is generally (1-4):1 in an indoor population (Kuang & Cheng, 1990). b The body color of the 17 females was an unusual wax yellow, which was different from the typical adult female color of T. cinnabarinus and T. urticae (the typical adult female color for T. cinnabarinus was brownish red and that for T. urticae was light green) ( Figure S2).
of adult female mites to these acaricides is shown in

| The effect of acaricides on mixed species
Mixed populations of the two mite species were generated in the laboratory for simulating the effect of acaricides on the population composition in the field. In the control group (water treatment), the ratio of Tu-YN to Tc-YN was 1 (5:5) at the beginning of the mixed population (origination) and reduced to 0.89 after 15 days; however, it dramatically declined to 0.46 two months following origination ( Table 6)

| The activities of detoxification enzymes (CarE, GST, and P450) in Tc-YN and Tu-YN
Compared with Tc-YN, Tu-YN has higher activities of detoxifying enzymes with the exception of GSTs. In Tu-YN, the P450 activity toward nitroanisole and the CarE activities were 1.81-and 2.80-fold higher than those of Tc-YN, respectively, whereas GST activity was not significantly different between the two species (Figure 3).    *The logarithm concentration-probit regression line passed the chi-square test. Chi-squared distribution was used to test whether the regression between logarithm concentration and probit pass the goodness of fit, if passed (the calculated chi-square values < expected χ 2 (df)0.05 , goodness-of-fit chisquare is significant), meaning the current data represent the fact.

| Induction effects of acaricides on activities of detoxifying enzymes in Tu-YN and Tc-YN
T A B L E 5 Susceptibility of Tetranychus cinnabarinus and Tetranychus urticae larvae to selected acaricides changes for P450s, GSTs, and CarEs in Tc-YN were 1.04-, 2.32-, and 0.90-fold, and 1.17-, 1.61-, and 1.07-fold in Tu-YN, respectively (Figure 4c). By and large, the responses of the three detoxification enzymes in Tu-YN were stronger than those in Tc-YN following the acaricide treatment; that is, the change in the activity of P450s in Tu-YN was significantly greater than it was in Tc-YN for all three treatment groups; furthermore, the increase in the activity of CarEs was higher for Tu-YN than Tc-YN for the abamectin-treated group, and the increase in the activity of GSTs was greater in Tu-YN than Tc-YN for the fenpropathrin-treated group, but lower for the tebufenpyrad-treated group, and where the changes in enzyme activity showed no significant differences between Tc-YN and Tu-YN in other circumstances (Figure 4a-c). filtering low-quality reads, the total numbers (mean) of clean reads of pre-exposure, abamectin-, fenpropathrin-, and tebufenpyradexposed mites were 13.7, 13.0, 13.1, and 13.8 million in Tc-YN, respectively (Table S1). The total numbers (mean) of clean reads were 18.2, 16.5, 17.7, and 15.6 million in Tu-YN, respectively (Table   S1). For each sample, over 97% of total clean reads were successfully mapped onto the T. urticae genome without mismatch for further analysis (Table S1). Tetranychus cinnabarinus and Tetranychus urticae are two closely related species, which facilitates us taking the advantage of the sequenced T. urticae genome and using it as a reference.

| The results of RNA-Seq for Tu-YN and Tc-YN
For interlibrary comparison, read numbers were normalized to relative abundance as reads per kilobase transcriptome per million mapped reads (RPKM). The RPKM value (mean value of three biological replicates) of each gene was further used to compute the related coefficients between each sample. The expression correlations of genes showed a better accordance between the same species under different acaricide exposures, whereas the correlations were markedly lower between Tc-YN and Tu-YN, both with and without exposure to an acaricide ( Figure S3). All raw and processed data have been deposited in the NCBI Gene Expression Omnibus under accession number GSE75529.  Figure S4).

| Functional annotation of differentially expressed genes (DEG) by GO, COG, and KEGG
All DEGs were aligned to the clusters of orthologous groups (COG) database for functional prediction and classification. A total of 2,300 COG annotations were identified for 253 annotated DEGs, which were classified into 24 molecular families ( Figure S5). Thus, some of these DEGs were associated with multiple COG annotations. Among these functional classes, the "general function prediction only" cluster constituted the largest group (131; 51.78%), followed by "carbohydrate transport and metabolism" (97; 38.34%), and "intracellular trafficking, secretion, and vesicular transport" (77; 30.44%).
To identify the metabolic pathways populated by these DEGs, all DEGs were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. A total of 587 DEGs were annotated to 198 KEGG pathways (p-value≤.05). The pathways with the most DEGs were "metabolic pathways" (116; 19.76%), followed by "lysosome" (47; 8.01%) and "retinol metabolism" (37; 6.3%). The pathway analysis showed that over 60% of DEGs are closely linked to the metabolism of xenobiotics and endogenous compounds (Table S2).

| Focus on DEGs involved in insecticide detoxification
The DEGs relating to insecticide resistance generally encode detoxification enzymes. Generally, P450s, CCEs, and GSTs are the three primary enzymes involved in the detoxification of insecticides.
Compared with Tc-YN, the gene expression profiles of P450 (CYPs) genes revealed that 15 CYP genes showed significant transcription level variations (10 genes upregulated and five genes downregulated) in Tu-YN ( Figure 5). These differentially expressed CYP genes were distributed among all CYP clans (clan 2, 3, 4, and M). The majority of the upregulated CYP genes (6 of 10) were represented by clan 2, including tetur03 g09961 and tetur23 g00260, which both had an over 10-fold increase in gene expression (Figure 9, Table S3). Three upregulated CYP genes belonged to clan 4, and another one was grouped to clan 3 (Table S3). comparing Tc-YN with Tu-YN, there were six upregulated and four downregulated GST genes, with three of the upregulated GST genes belonging to class mu and three of these genes belonging to class delta ( Figure 5, Table S3). In this study, we found that 13 and 7 CCE genes were significantly up-and downregulated, respectively, in Tu-YN vs. Tc-YN ( Figure 5, Table S3).
There were several other biotransformation enzymes potentially participating in the detoxification process, such as short-chain dehydrogenase/reductase (SDR), sialin, ABC transporters, and intradiol ring-cleavage dioxygenase (ID-RCDs). The genes that encode for these enzymes were classified as "others" in this study to separate them from the P450s, GSTs, and CCEs. In Tu-YN vs. Tc-YN, 10 "others" genes (five SDR genes, one sialin gene, three ABC-transporter genes and one ID-RCDs gene) were upregulated, and four (three SDR genes and one  Table S3). In total, there were 39 and 20 genes involved in biodegradation or transport (including P450s, GSTs, CCEs, and "others") that were up-and downregulated, respectively, in Tu-YN compared with Tc-YN. This showed that the number of specifically upregulated detoxification and transport genes in Tu-YN was nearly twice that of Tc-YN (Figures 5 and 9).

| Comparative analysis of the response to acaricides between Tc-YN and Tu-YN using RNA-Seq
To obtain a global view of transcriptome responses following treatment with acaricides, we analyzed variations in the gene expression profile between acaricide-treated and untreated mites (as a negative control group) using digital gene expression (DGE), which is a highthroughput tag sequencing (Tag-seq) method used to identify up-and downregulated genes between two datasets. To identify differentially expressed genes following exposure to acaricides, differences in gene expression were analyzed by pairwise comparisons between mites with and without acaricide exposure. We found that 35 and 133 genes were significantly up-and downregulated, respectively, in the control-vs. abamectin-exposed in Tc-YN, while in Tu-YN 89 and 179 genes were up-and downregulated following abamectin exposure (Figure 6a). Among these differentially expressed genes, 73 genes were coregulated (22 upregulated and 51 downregulated) in both species following abamectin exposure. Among the differentially expressed genes in the two comparison groups, one P450 gene, three GSTs, and three "others" genes were upregulated in both Tc-YN and Tu-YN; however, seven P450s, six GSTs, two CCEs, and six "others" genes were upregulated in Tu-YN, whereas only one P450, five GSTs, one CCE, and three "others" genes were upregulated in Tc-YN.
Following fenpropathrin exposure in Tc-YN, 32 and 33 genes were significantly up-and downregulated in the treatment group compared with the unexposed group, while in the exposed Tu-YN group there were 74 upregulated and 104 downregulated genes compared with the untreated Tu-YN controls (Figure 7a). Among the differentially expressed genes in the two comparison groups, 11 P450s and five "others" genes were specifically upregulated in Tu-YN (Figures 7b and   9), while only one P450 and one "others" gene were upregulated in Tc-YN. Of the remaining detoxification enzyme categories, two GSTs genes and two CCE genes were specifically upregulated in Tu-YN and Tc-YN, respectively (Figure 7b, Table S5).
Following tebufenpyrad exposure, there were 21 up-and 50 downregulated in Tc-YN, whereas 59 genes were upregulated and 120 were downregulated in Tu-YN (Figure 8a, Table S6). Of the upregulated genes in Tc-YN, none were P450s, GSTs, or CCEs; however, three P450 and four GST genes were found to be upregulated in Tu-YN (Figures 8b   and 9). Furthermore, six "others" genes were specifically upregulated in tebufenpyrad-treated Tu-YN, while only one gene within the "others" category was upregulated in Tc-YN (Figures 8b and 9).
Notably, a specific enrichment of genes belonging to certain metabolism pathways was observed in all DEGs, such as "metabolic pathways," "metabolism of xenobiotics by cytochrome P450," "drug metabolism-cytochrome P450," "retinol metabolism," "arachidonic acid metabolism," "ascorbate and aldarate metabolism," and "drug metabolism-other enzymes" (Tables S7-S9). The other notable phenomenon was that the DEGs that belonged to certain metabolism and sequestration pathways in Tu-YN were much more abundant than those in Tc-YN following acaricide exposure (Tables S7-S9).

| Sequence alignment of target of acaricide genes in Tu-YN and Tc-YN
A total of 36 SNPs (nine nonsynonymous) were identified for the 14 acaricide-target genes in Tu-YN, while 32 SNPs (six nonsynonymous) were observed in Tc-YN (Table 7). While multiple SNPs were identified for both mite species, none of the SNPs were located within the active site of the proteins encoded by these genes nor reported as a resistance-related in previous studies (Ilias et al., 2014;Kwon et al., 2015). These results suggested that the identified SNPs are not correlated with acaricide resistance.
F I G U R E 6 Summary of differentially expressed genes (DEGs) in each pairwise comparison following abamectin exposure. (a) Venn diagrams depicting overlap among DEGs of Tc-YN and Tu-YN exposed to abamectin, respectively. For each Venn diagram section, the numbers of transcripts differentially expressed in any strain with treatment as compared to control were indicated. (b) Distinctive and common differentially expressed (upregulated) genes of detoxification enzymes in the two previous comparisons. Tc-AV: Tc-YN were treated with abamectin; Tu-AV: Tu-YN were treated with abamectin. The category "others" included SDR, ID-RCD, and ABCtransporter genes

| DISCUSSION
The differentiation between T. urticae and T. cinnabarinus is often difficult and controversial; however, in China, this is simplified in practice. That is, the green (female adult) with two black brown spots is named the two-spotted spider mite and viewed as an invasive species, while the red (female adult) with two black brown spots is named the carmine spider mite and viewed as a native species (Kuang & Cheng, 1990;Sun et al., 2012). For more accuracy, the number of setae on tibia I and cross with T. cinnabarinus maintaining in laboratory were studied for identifying the species of collected field population in this work.
According to Kuang and Cheng (Kuang & Cheng, 1990), the number of setae on tibia I for T. urticae population is held constant at 10, while in T. cinnabarinus the majority of individuals have 10 with a smaller number having 12 or 13. In this study, the proportion of individual with 10 setae on tibia I was 100% in the T. urticae field population, while it was only 70 and 66.7% in the T. cinnabarinus field and laboratory populations, respectively. In addition, the cross between the T. cinnabarinus and laboratory population could continuously produce F 1 , F 2 , and could bias the results. However, the method is adequate to confirm that T. urticae has expanded successfully, as it has been reported from more collection locations than T. cinnibarinus, despite having originally invaded only in Beijing. We suggest that T. cinnabarinus is contracting its distribution based on two facts. First, the practice of in-field observations from ourselves and peers found that the collection of T. urticae is presently much easier than that for T. cinnabarinus; second, F I G U R E 9 Clustering of metabolic enzyme genes differentially expressed across strains. Hierarchical clustering analysis-based transcription levels was performed on 81 enzyme-encoding genes showing significant differential transcription (Fisher's test p value <.001) in Tu-YN compared with Tc-YN for samples both exposed to acaricides and compared to their respective controls. Clustering was obtained using Euclidean distance calculated from all Log 2 (fold changes) and complete linkage algorithm. Color scale from green to red indicates Log 2 (fold changes) from -2 (fourfold under transcription) to 2 (fourfold over transcription). For each gene, either the gene ID from the T. urticae genome or the functional annotation of the gene is indicated. (a) Tu-YN vs. Tc-YN, (b) Tc-YN exposed to abamectin vs. Tc-YN, (c) Tu-YN exposed to abamectin vs. Tu-YN, (d) Tc-YN exposed to fenpropathrin vs. Tc-YN, (e) Tu-YN exposed to fenpropathrin vs. Tu-YN, (f) Tc-YN exposed to tebufenpyrad vs. Tc-YN, (g) Tu-YN exposed to tebufenpyrad vs. Tu-YN. tetur06 g04880 tetur106 g00040 tetur06 g05080 tetur06 g05090 tetur28 g01720 tetur02 g01350 tetur28 g01130 tetur12 g02130 tetur07 g07310 tetur22 g00170 represent RCD gene; tetur01 g00490 tetur06 g00460 tetur04 g08620 tetur28 g01250 represent ID-RCD; tetur02 g09080 represent sialin; tetur07 g04290 tetur07 g04410 tetur03 g09880 tetur06 g03560 tetur06 g03510 represent ABC transporters. FC: Fold Change a 5-year-long survey has demonstrated that T. urticae has supplanted T. cinnabarins as the dominant species of mites on vegetables in Beijing and Hebei provinces (Wang et al., 2013). All of the above-provided indirect evidence that the range of the invasive T. urticae is expanding and outcompeting the native T. cinnabarinus in China.
In ecological theory, when one species show a competitive advantage against others in the same environmental conditions, the advantage could be attributable to the better ability of the more advantageous species to fit to the environment (Reitz & Trumble, 2002;Roush & Daly, 1990), also known as a fitness advantage.  (Bolland et al., 1998;Saito, 1979); however, the competitive capacity was found to vary between T. urticae and T. cinnabarinus when they fed on different host plants (Li & Cheng, 2011;Liu & Sun, 1998;Tomczyk, Kropczynska, & Elst, 1995). Theoretically, the distribution of host plants may be an extremely important factor in which one species has the overall advantage; nevertheless, we could not judge the role of host plant type on the distribution of these two species due to the lack of relevant data from investigation to date. In addition, temperature may be excluded as an impact factor as T. urticae is typically distributed in high latitude regions, which are characterized by lower temperatures, while T. cinnabarinus prefers regions of lower latitude (Dupont, 1979;Goka & Takafuji, 1991;Takafuji, So, & Tsuno, 1991). The higher fecundity of T. urticae at lower temperatures (16°C) and lower fecundity at higher temperatures (33°C) compared with T. cinnabarinus were also confirmed in our study (Table 3). We can imagine from their nature that T. cinnabarinus would benefit from global warming and enlarge its distribution wider than T. urticae; however, this is not true (Figure 2c). Over the past three decades, T. urticae gradually spread into the south of China; that is, T. urticae spread from high latitude regions of lower temperatures to low latitude regions of higher temperatures in China. Therefore, the factor of temperature cannot be considered crucial for the competitive capacity of T. urticae expansion to become the dominant spider mite in China.
Existing studies conformably revealed that T. urticae is less sensitive to most acaricides than T. cinnabarinus and that the difference in sensitivity to abamectin, an acaricide with extensive use in China, is especially significant (Bi et al., 2016;Gu et al., 2000;Zhao et al., 2006).
Similar bioassay results were obtained in our study (Tables 4 and 5).
These highly consistent results encourage us to make a bold inference -32 6 a VGSC: voltage-gated sodium channel gene; GluCl-01~GluCl-05: 5 glutamate-gated chloride channel genes; Rdl1-Rdl3: 3 GABA receptor genes; SdhA-SdhE: 5 succinate dehydrogenase complex genes. b Alignment between tested cDNA sequence and baseline at the level of nucleotide. c Alignment between tested cDNA sequence and baseline at the level of amino acid. d Identified point mutation from amino acid sequence was matched with the reported mutations of target genes involved in acaricide resistance in mites; "-" indicates not match.
that it is the higher tolerance against acaricides, in particular against abamectin, which is the major factor that has resulted in the expan-  (Chu, Wan, Zhang, & Brown, 2010;Guo et al., 2014).
This difference allows the Q biotype to overcome the competitive advantage that the B biotype has in the absence of insecticide pressures (Elbaz et al., 2012;Pascual & Callejas, 2004). Similarly, Gao, Lei, Abe, and Reitz (2011), Gao et al. (2012) suggested that differential susceptibility to commonly used insecticides could account for the replacement of Liriomyza sativae by Liriomyza trifolii on Hainan Island of southern China (Gao et al., 2011(Gao et al., , 2012. However, the underpinning mechanisms behind the differential susceptibility to insecticides for these Bemisia spp. and Liriomyza spp. were not studied in depth or elucidated. It is reasonable for us to compare the activities of P450s, GSTs, and CCEs between T. urticae and T. cinnabarinus as they are the main enzymes functioning in the detoxification and metabolizing of exogenous chemicals, such as a variety of insecticides (Bass & Field, 2011;Li, Schuler, & Berenbaum, 2007;Qiu, 2005). Our results showed that the activity of GSTs did not differ significantly between T. urticae and T. cinnabarinus; however, the activities of P450s and CarEs were significantly higher (1.81-and 2.80-fold, respectively) in T. urticae than in T. cinnabarinus. Furthermore, the change in enzyme activity showed a significant difference between T. cinnabarinus and T. urticae in six of nine cases of in-laboratory acaricide-stimulated experiments, in which T. urticae responded more strongly in five cases while T. cinnabarinus responded more positively in only one case (Figure 4). These results from the enzyme activity studies revealed from a biochemical perspective that T. urticae possesses a stronger capacity for the detoxification of acaricides than T. cinnabarinus.
At a molecular level, the transcriptome contains the complete repertoire of mRNAs transcribed by a living cell, that is, the sum of genetic information transcribed from the genomic DNA .
Comparative expression profiling between T. cinnabarinus and T. urticae provides a priori information that can reveal the underlying mechanisms mediating the stronger acaricide tolerance of T. urticae over  (Tables   S7-S9). All of the above results strongly suggested that compared with T. cinnabarinus, T. urticae possesses a greater potential and ability to mitigate the stress from acaricides by possessing more upregulated genes (including metabolism-related genes) at a basal gene expression level, and more genes that respond by intensively increasing their expression levels following acaricide exposure. This is the underlying molecular reason that T. urticae possess a stronger tolerance to acaricides than T. cinnabarinus. An important thing need to point out is the detoxifying genes switching them on more in the presence of three acaricides compared with no acaricide exposure, which may be explained that higher expression of detoxifying genes would come at a cost and it may only be beneficial to express these genes more when a pesticide is present.
In addition to our investigation of well-known detoxification gene families involved in acaricide resistance, our study on target-site modification revealed that the number of identified SNP differences within the acaricide-target genes from the two field populations was not different too much. Furthermore, none of the identified SNPs were consistent with the reported mutations of target genes known to confer target-site resistance against pesticides (Ilias et al., 2014;Kwon et al., 2015). The sequence alignment of target genes in T. urticae and in T. cinnabarinus implied that target-site modification probably played a very minor role in conferring the higher acaricide tolerance in T. urticae over T. cinnabarinus.
In summary, this paper provides evidence that the activities and expression levels of detoxification enzymes were generally greater in T. urticae than in T. cinnabarinus. This enhanced metabolic detoxification may be an important reason why T. urticae is more resistant than T. cinnabarinus to acaricides. Furthermore, this greater resistance to acaricides may explain the phenomenon that acaricides facilitate the continuous expansion of T. urticae as the dominant spider mite in many locations in China, excluding other aspects as an essential factor. The population competition experiments further supported this hypothesis that the competitive displacement of T. cinnabarinus by T. urticae is mediated by human activities. Not only do our results reveal that T. urticae possesses stronger detoxification capacity than its sibling species T. cinnabarinus, which facilitated its persistent expansion in China, but they also point to the need to accurately identify Tetranychus species and to develop species-specific management strategies for these pests.