Species delimitation of the Dermacentor ticks based on phylogenetic clustering and niche modeling

Sample collection

D. silvarum and Haemaphysalis longicornis were collected from Xiaowutai National Natural Reserve of Hebei Province, China. D. nuttalli were collected from a farm in the suburb of Yanbian Korean Autonomous Prefecture, Jilin Province, China. All ticks were first identified morphologically under a light microscope, and then verified by molecular data.

DNA extraction, PCR and sequence analysis

Total genomic DNA was extracted from ticks using the DNeasy blood and tissue Kit (Bioteke, Beijing, China) following the manufacturer’s protocols. Amplification of COI and ITS2 were performed using the respective primer pairs: LCOI490F (5′- GGTCAACAAATCATAAAGATATTGG-3′)/HCO2198R(5′-TAAACTTCAGGGTGA- CCAAAAAATCA-3′) (Folmer et al., 1994; and F3/1 (5′-GGGTCGATGAAGAACGCA- GCCAGC-3′)/R1/1 (5′-TTCAGGGGGTTGTCTCGCCTGATG-3′) (Kulakova et al., 2014). The PCR reactions contained a total of 40 µl consisting of 20 µl 2 ×Taq PCR Master Mix (Biomed, Beijing, China), 17 µl distilled water, 1 µl of each primer and 1 µl DNA template, and a negative control was included for all reactions. The amplifications were conducted on the ProFlex™ PCR System (Thermo Fisher, Waltham, MA, USA). The cycling conditions for COI was as follows: 94 °C for 5 min, 35 cycles of 94 °C for 30 s, 43 °C for 50 s, 68 °C for 1 min, and an extension of 68 °C for 10 min. The protocols for ITS2 was: 95 °C for 10 min, 30 cycles of 94 °C for 30 s, 58 °C for 30 s, 72 °C for 1 min, and an extension of 72° C for 10 min. The PCR products were visualized on 1% agarose, and the most intense products were sent for sequencing in Sangon Biotech (Shanghai, China). Moreover, sequences of COI and ITS2 of portion ticks were obtained from GenBank (http://www.ncbi.nlm.nih.gov/genbank/).

All the sequences were analyzed using Mega 7 (Kumar, Stecher & Tamura, 2016) and phylogenetic analyses were conducted using MrBayes v3.1.2 (Ronquist & Huelsenbeck, 2003). The HKY+I and GTR model were respectively selected for COI and ITS2 data using MrModelTest 2.3 (Nylander, 2004). For the Bayesian analyses, two independent runs were performed for 5,000,000 generations by sampling one tree per 1000 generations. The first 25% of trees were discarded from analysis as burn-in. H. longicornis (Ixodidae) was chosen as the outgroup.

Digital occurrence records

The distribution information of ticks was collated from the Global Biodiversity Information Facility (GBIF) database (https://www.gbif.org/) and Vectormap (http://vectormap.si.edu/index.htm). Further data was obtained from the related published literature (Table S1) and from the georeferenced specimens in the Hebei key laboratory of Animal Physiology, Biochemistry and Molecular Biology, China. Online gazetteers and Google Earth were used to reconstruct the geo-coordinates for each chosen datapoint. Occurrence records of each species were double-checked by DIVA-GIS software (version 5.2) to detect possible errors in georeferencing. A total of 685 occurrence records for Dermacentor species were assembled in our analysis (Table S2). Occurrence records are often biased due to the different densities of species distribution. Thus, in order to reduce the sampling bias and remove the spatial autocorrelation, we used a coarse resolution of 5 arc-minute and then randomly selected a single point in each grid cell (Rodríguez-Castaneda et al., 2012; Li, Du & Guo, 2015).

Climate variables

To characterize the climate heterogeneity across the distribution range for the tick species, we initially compiled all 19 bio-climatic variables, which were download from the Worldclim Global Climate Database (http://www.worldclim.org) (Hijmans et al., 2005), representing minima, maxima and average values of monthly, quarterly, and annual ambient temperature as well as precipitation recorded between 1950 and 2000. Moreover, a spatial resolution of 2.5 arc min (approx. ∼5 km resolution at the equator) of these climate data was used.

The accuracy of a model can be affected by the strong covariance among environmental variables. In order to minimize the multicollinearity among predictor variables, we used principal component analysis (PCA) of all 19 climatic variables in SPSS statistics to identify and remove highly correlated variables (—r— > 0.80) from our models (Table S3). Finally, six climatic environmental variables were selected for the tick species: Bio2 (Mean Diurnal Range), Bio8 (Mean Temperature of Wettest Quarter), Bio12 (Annual Precipitation), Bio15 (Precipitation Seasonality), Bio18 (Precipitation of Warmest Quarter) and Bio19 (Precipitation of Coldest Quarter).

Climatic niche modeling

In the last few decades, species distribution models (SDMs) became a popular tool for estimating the potential distribution for animals (Park et al., 2017; Wei et al., 2018) and plants (Thuiller et al., 2005; Yi et al., 2016). Lots of SDMs were designed to predict potential distribution of a species, such as MaxEnt, GARP and BIOCLIM (Beaumont, Hughes & Poulsen, 2005). However, the MaxEnt software (version 3.3.3k, http://www.cs.princeton.edu/sc˜hapire/maxent/) has been used mostly frequently to simulate species potential distribution, since it only requires presence-only data (Radosavljevic & Anderson, 2013; Cooper et al., 2016). Many studies had shown that MaxEnt performed well regardless of the number or geographical extent of species records, and outperformed than other methods when predicting the potential distribution of species, especially for the small individual samples (Hernandez et al., 2006; Elith et al., 2011; Bosso et al., 2016; Sultana et al., 2017).

Some recent studies had demonstrated the importance of MaxEnt model settings to balance the performance and complexity of the model, as the default settings of MaxEnt could result in overfit models (Radosavljevic & Anderson, 2013; Merow et al., 2014; Warren et al., 2014). Thus, to avoid overfitting, simultaneously maximize the predictive power and provide the best model for the species, we adopted the R package (version 3.4.2, R Core Team, 2017) ENMeval to select the optimal combination of two important MaxEnt parameters, the value of the regularization multiplier and the combination of feature classes (Muscarella et al., 2014; Naveda-Rodríguez et al., 2016). The “checkerboard2” approach was employed to calculate the standardized Akaike Information Criterion coefficient (AICc). The parameterizations that resulted in the model with the lowest delta AICc score were selected to run the final MaxEnt models (Warren & Seifert, 2011; Radosavljevic & Anderson, 2013). As AICc can only choose the ‘best’ one from a set of models, and can not directly assess the model performance, we inspected the omission rate and tested AUC (area under the curve) of the models that were selected as optimal. Moreover, the parameter for regularization multiplier and feature were set based on the ENMeval analysis. The regularization multiplier was varied from 0.5 to 4 in increments of 0.5, and the following five feature classes were tested: (1) Linear (L); (2) Linear (L) and Quadratic (Q); (3) Linear (L), Quadratic (Q) and Hinge (H); (4) Linear (L), Quadratic (Q), Hinge (H) and Product (P); (5) Linear (L), Quadratic(Q), Hinge (H), Product (P) and Threshold (T).

A minimum convex polygon that was comprising all records of each species were randomly chosen to define background points (Phillips, 2008). Additionally, the logistic output of MaxEnt was used for all analyses. The LPT (lowest presence threshold) was used to define the suitable and unsuitable habitats for all three species. This threshold was a conservative value that was widely used in species distribution modeling, especially when data was collected by different observers and methods over a long period of time (Bosso et al., 2016). The final model was converted into a binary presence/absence map, which was created using the reclassify module from ArcGIS 10.2 (ESRI) by applying this threshold. To maximize the predictive information and simplification of future analysis, the suitable habitat areas for the three tick species were classified into four levels.

Modeling evaluation

The AUC of the receiver operating characteristic (ROC) (Fielding & Bell, 1997) were used to estimate the performance of the model. The AUC value was ranged from 0 to 1, where a value below 0.5 was interpreted as a random prediction; 0.5–0.7 indicated poor model performance; 0.7–0.9 indicated moderate performance; and a value above 0.9 was considered to have ‘good’ discrimination abilities (Peterson et al., 2011). A 10-fold cross-validation was used to run MaxEnt to prevent random errors from affecting the selection of the validation and prediction samples. To assess the influence of environmental variables on species, the jackknife test was adopted to measure the importance and the percent contributions of each variable.

Niche overlap test

Niche overlaps among species were quantified using the PCA-env method proposed by Broennimann et al. (2012). PCA was used to transform the environmental space of the investigative environmental variables into a two-dimensional space defined by the first two principal components (Strubbe et al., 2013; Strubbe, Beauchard & Matthysen, 2015). Then, the two-dimensional environmental space was projected onto a 100 ×100 PCA grid of cells bounded by the maximum and minimum PCA values in the background data. A kernel function was used to smooth the climatic space defined in the gridded PCA climatic spaces based on the first two principal components (PCs) (Petitpierre et al., 2012). The niche overlap among the species was measured directly by Schoener’s D from the ecological niche space (Warren, Glor & Turelli, 2008). The Schoener’s D was an index which varied from 0 (no overlap) to 1 (overlap). Additionally, niche equivalency and a similarity test were also implemented and if the niche overlap value fell outside the 95% confidence interval of the null hypotheses, equivalency of the two niches was rejected. And in the similarity test, a p value >0.05 indicated that the niches were no more similar than expected by chance.

All GIS analyses were performed using ArcGIS version 10.2 (ESRI). The “ecospat” package in R was used to implement this analysis (Broennimann, Di Cola & Guisan, 2018).

Phylogenetic divergence

Nucleotide sequences obtained in this study were used to assess the phylogenetic relationship of the tick species. The length of the COI sequences were 625 bp after edge trimming. ITS2 sequences were ranged from 857 bp in H. longicornis to 1039 bp in Dermacentor species. Two phylogenetic trees (based on COI and ITS2) obtained from Bayesian analyses showed similar topology (Figs. 1A, 1B). All obtained sequences clustered tightly in a branch and formed a distinct lineage. Moreover, the two species D. silvarum and D. nuttalli grouped together, and D. marginatus was in a separate cluster in the two analyses.

ENMeval data and model performance for potential distribution

The results of ENMeval were shown in Table S4 and Fig. S1. Regularization multiplier = 1, feature combinations = Linear, Quadratic, Hinge, Product and Threshold (LQHPT) were chosen to MaxEnt configuration for D. marginatus; regularization multiplier = 0.5, feature combinations = Linear and Quadratic (LQ) were chosen as the optimal configuration for D. nuttalli and D. silvarum.

Based on 10-fold cross validation, the model performance for the three species were all better than random, with a mean AUC value of 0.948 ±  0.013 for D. marginatus, 0.952 ±  0.019 for D. nuttalli and 0.940 ±  0.044 for D. silvarum. These results indicated that the model performed well in predicting the species distribution for all three species. A threshold value of 0.0021 for D. marginatus, 0.0034 for D. nuttalli and 0.0037 for D. silvarum, which were obtained from the LPT analysis. The suitable habitat area for these species were classified into four levels: <threshold indicated unsuitable habitat; threshold-0.4 indicated low habitat suitability; 0.4–0.6 indicated moderate habitat suitability; and >0.6 indicated high habitat suitability.

Important environmental variables for ticks

The relative contributions of the six environmental variables for the three species were shown in Table 1 and Fig. 2. Precipitation of coldest quarter (Bio19) and Mean temperature of wettest quarter (Bio8) had the largest contributions to the distribution model for D. marginatus. These two factors explained 73.3% of the modeled distribution. The distribution of D. nuttalli was mainly constrained by precipitation seasonality (Bio15) and Mean temperature of wettest quarter (Bio8), which explained 54.6% of the distribution. D. silvarum was mainly constrained by Precipitation of warmest quarter (Bio18) and Mean temperature of wettest quarter (Bio8), which accounted for 53% of the modeled distribution. Mean temperature of wettest quarter (Bio8) was the second highest contribution for the model distribution in all three species. Among the six environmental variables that our analysis focused on, precipitation conditions were more important than other factors when creating the distribution map for these species. The potential distribution of these species were shown in Fig. 3 and Fig. S2. The potential distribution maps predict that D. marginatus would be found mainly in European, such as France, Italy, Spain, Greece and Portugal, which have highly suitable climate for this species. However, our analysis predicted that D. nuttalli and D. silvarum would have similar distribution area which mainly focus on eastern of Eurasia. Moreover, compared to D. silvarum, D. nuttalli had wider potential distribution, including Mongolia and northern China.

Niche comparisons

The ordination approach using PCA-env revealed the niche patterns for each species pair in the environment space (Table 2 and Fig. S3). The PCA-env resulted in the species pair D. marginatus and D. nuttalli showed two PCA axes explaining 71.98% of all climatic variables, and two PCA axes explaining 73.24% for species pair D. marginatus and D. silvarum. The analysis also revealed that 74.28% of all climatic variables for D. nuttalli and D. silvarum were explained by two PCA axes. Furthermore, niche overlaps showed a great variability in the E-space inhabited by the different species. There were 17% niche overlap between D. nuttalli and D. silvarum (Schoener’s D = 0.17), 9% between D. marginatus and D. nuttalli (Schoener’s D = 0.091), and only 2% between D. marginatus and D. silvarum (Schoener’s D = 0.026). Of the three species that were analyzed, niche overlaps of D. nuttalli and D. silvarum were highest.

The environmental space occupied by each species as determined by PCA-env were also shown in Fig. S3 and Table 2. The null hypothesis of the niche equivalency test was rejected for all possible pairwise comparisons, suggesting that climate niche among these species pairs were significantly distinct. On the other hand, in niche similarity analysis, null hypothesis testing were held on all species pairs. For the Dermacentor species investigated, the niche similarities were higher than expected by chance. Taken together, we found that the niche of these three species were similar but not identical, highlighting that while these three species are closely related, they represent distinct species.