Species distribution modeling
Species data were simple presence observations based on records from checklist and our field sampling. Occurrence records for D. dubium were mostly referred to [22], our consecutive sampling as well as field observation. Occurrence records for D. excisum were mostly from the published literature. D. excisum was widely distributed in tropics and subtropics of Australia, Asia and Africa, up to eastern China. No records were reported for both species from South America, North America and Antarctic. In total, 165 populations for D. dubium and 162 populations for D. excisum were used for species distribution models (see Additional file 1).
We used a comprehensive set of 19 bioclimatic variables with a spatial resolution of 2.5 arc min (about 4 km in the study area) as climate predictors available through http://www.worldclim.org [26]. We computed pair-wise spearman rank correlations to assess multi-colinearity of the 19 wordclim variables, wherein species distribution models (SDMs) were developed with a subset of variables with R2< 0.75: bio2, bio3, bio5, bio6, bio8, bio13, bio14, bio15, bio18, bio19 (see Additional file 1). For each species we use the modelled association between current climate and present-day distribution to estimate current distributional areas.
We performed spatial predictions of the two species distributions in ‘ecospat’ written in R language [27], which provided methods and utilities for quantifying climatic niches shifts between two congeneric peripheral distribution species, analyzing niche overlap between them [28]. We used a shape file containing presence-absence records for a species as spatial points generated in QGIS v2.18. Raster datasets (in grid format) of above ten wordclim variables were used as explanatory variables (predictors). Coupled with selected explanatory variables, the shape files with only occurrence data of both D. excisum and D. dubium were used to perform modeling and evaluation with package ‘sdm’. Four models (‘rf’ (Random Forest), ‘brt’ (Generalized Boosting Model or usually called Boosted Regression Trees (gbm)), ‘bioclim’ (Surface Range Envelop (sre), ‘glm’ (Generalized Linear Model)) were used for model training. Roc curves were used to compare the results for all models. Then we just predicted the habitat suitability into the whole study area, and fitted habitat suitability value of each present-absent site for each model were obtained from the output shapefiles in QGIS v2.18 (see Additional file 1).
For pre-modeling, we investigated spatial autocorrelation of climate predictors by Mantel correlogram to determine the minimal distance between species records. We performed Niche quantification and comparison with ordination techniques. The PCA was calibrated on all the sites of the study area to describe the contribution of all original predictors, covering the ranges of the two study species [29]. Niche overlap in terms of Schoener’s D was computed based on kernel density estimates for both species occurrence records in PCA space accounting at the same time for the relative frequency of specific climate conditions within the study area. The results allowed a quantification of those parts of niche space occupied by both species as well as the unique parts of two species. The classification of Schoener’s D was proposed by Rödder and Engler [30].
Randomization tests for niche equivalency and similarity as proposed by Warren et al. [31] were computed to test the hypothesis of niche equivalency and niche similarity [32]. Parameter alternative was set as ‘greater’ to test for niche conservatism, i.e. whether the overlap is more equivalent/similar than random. For Core Niche Modeling, we accessed the calibration capacity of our models using a repeated split-plot approach [33]. We used a random subsample of 75% of the points of presence to calibrate the models using the approach described above, and the remnant 25% to evaluate the model, repeating the whole procedure 10 times. Presence-only data were used to calculate the Boyce Index for each model [34], as Boyce index measured how much model predictions differ from a random distribution of the observed presences across the prediction gradients [35]. Finally, we ensemble selected small models with the Boyce index more than 0.80 for spatial predictions with current climatic variables. Once the models were calibrated and evaluated, we projected the potential distribution of the species over space. Spatial projections were visualized and plotted in QGIS v2.18. All statistical analyses were performed using R 3.5.3 [27], and the packages being “biomod2”, “ecospat”, “sdm’, “maptools”, “raster”, “rgdal”, “rtiff”.
Life history and temperature niche
Analysis of temperature niche difference between the two species was conducted with two representative clones, one for D. dubium and another for D. excisum. D. dubium (clone N) was collected from a fish pond on the way to Changtang reservoir, Shaoguan, China (N24°78’, E113°50′). D. excisum (clone D) was collected from Hedi reservoir, Zhanjiang, China (N21°85’, E110°31’). The two clones chosen were each the best performing clone of its species from our previous life history experiments with 16 Diaphanosoma clones [25]. We measured their life history traits (net reproductive rate (R0), intrinsic rate of population increase (r)) under a temperature range from 10°C to 40°C. 40℃ was estimated to be above the critical thermal maximum (CTmax), while 10℃was estimated to be below the critical thermal minimum (CTmin) [25]. Life table experiments were conducted at 15, 20, 25, 30, 34, and 38 °C. There were at least five replicates per clone by temperature combination; replicates being single individuals in 30 mL of medium. Additional replicates were added in environments with high juvenile mortality (T > 34℃). The detailed life history experiments and parameter calculation were described in Pajk et al. [25]. As temperature curve performance of R0 and r can be used as representative of temperature niche, they were the fitted by mixed linear models with the lme function from the R package “nlme” [36]. The choice of the best fit was made through comparison of the AICc.
Interspecific competition experiments
The competition experiments were conducted with above clones of two species of Diaphanosoma. The 12±12 h old second clutch offspring of the F3-generation were used to start the experiments. The experiments were run at two temperatures (20°C and 30°C) and food levels (1.0 mgC/L and 0.5 mgC/L) [37]. The temperatures chosen are representative of the localities of D. dubium and D. excisum, while the food levels are selected from the life history experiments of D. dubium and D. excisum, in which a sound intrinsic rate of increase (r>0.2) was obtained at above 0.5 mgC/L.
Animals were kept in 300 mL of medium. Monoculture replicates initially contained 10 neonates of a single species. Competition replicates were started with 10 neonates of each species. The number of animals in each treatment was recorded every six days at 20°C and every four days at 30°C. For each replicate we determined the maximum density and the maximum growth rate
N0 and Nt: population densities at two consecutive samplings). We calculated r of the two species in both competition cultures and monocultures [38-40]. These parameters were compared between the treatment groups (species×T×food×monoculture/competition culture) with a factorial ANOVA (SPSS 16.0).