Correlation of biotic and abiotic factors to insect community composition
We first used Mantel tests to assess the overall relationship between insect and tree (species and phylogenetic) β-diversity for plot and transect level matrices. Because Mantel tests perform poorly when accounting for the influence of geographic separation on patterns of association (Legendre et al. , 2015), we developed an ordination (redundancy analysis-RDA) approach to test tree and tree’s phylogeny combined with environmental variation and spatial distance to explain beetle composition. We used forward and backward selection in an RDA analysis assessing the influence of plant species (Hellinger transformed) and plant phylogenetic turnover on beetle composition (Hellinger transformed). The selection procedure was conducted using the ‘ordistep’ function in ‘vegan’. Variation partitioning analyses through redundancy analysis was used to assess the percentage contribution (both unique and shared) of each group of predictor variables (i.e., plant species composition/plant phylogenetic turnover, environmental, and spatial distance) to explain the variation in abundance-based longhorn beetle species composition. The environmental variables were converted to a distance matrix after log-normalization, while the spatial distance which was captured as latitude/longitude coordination was converted into the Cartesian coordination for the above calculation. Significance of testable fractions (P ≤ 0.05) were based on 999 permutations. Variation partitioning analyses were performed in R using ‘vegan’ and ‘packfor’ packages.
To further explore the potential mechanism of beetles community assembly pattern, nonmetric multidimensional scaling (NMDS) was conducted to show how dissimilarity changed within and between each sampling transect. NMDS analysis was performed using ‘metaMDS’ in ‘vegan’. The multi response permutation procedure and mean dissimilarity matrix (MRPP) algorithm was employed to ascertain if the mean distance within each group, which was defined as mean elevation of the transect, was significantly different from the mean distance of all the plots based on 999 permutations. To quantify the homogeneity of dissimilarity variances within each transect, the variances of the dissimilarity matrix were compared using the ‘betadisper’ method (Anderson, 2006). This test is analogous to Levene’s test for homogeneity of ANOVA variances. Finally, we used analysis of similarities (ANOSIM) to test if there was a significant difference between transects (Clark et al. , 1999; Warton et al. , 2012). The P -value was also obtained based on 999 permutations.
Finally, we introduced linear mixed-effect model to analyze the effect of plant α-diversity and phylogenetic α-diversity, environmental variability and geographical variability on beetle α-diversity, respectively. In total we considered 4 groups of datasets: set 1) beetles standardized Shannon diversity (BeeShannon_stdz); set 2) plant standardized α-diversity metrics, which including standardized Chao1 diversity (PlaChao1_stdz), standardized Simpson diversity (PlaSimpson_stdz), standardized Shannon diversity (PlaShannon_stdz) and standardized phylogenetic α-diversity (PlaPD_stdz); set 3) standardized environmental variability, which including standardized AMT (AMT_stdz), standardized ATR (ATR_stdz) , standardized AMH (AMH_stdz), standardized AHR (AHR_stdz), standardized MTWM (MTWM_stdz) and standardized MTCM (MTCM_stdz); set 4) standardized geographical variability, which including standardized ELE (ELE_stdz) and standardized Latitude (Latitude_stdz). Beetles standardized Shannon diversity were treated as response variable and site names were treated as random effect, and the remaining variables including set 2), set 3) and set 4) were treated as fixed effect. Moran’s I correlograms were built to evaluate the degree of spatial autocorrelation of the variables in relation to geographic distances and we found no significant positive spatial autocorrelation for these variables. For each dataset, we first fitted one global model including all the fixed effects and either a random intercept and slope (BeeShannon_stdz ~ PlaPD_stdz + PlaChao1_stdz +PlaSimpson_stdz + PlaShannon_stdz + ELE_stdz + AMT_stdz + ATR_stdz + AMH_stdz+ AHR_stdz + MTWM_stdz + MTCM_stdz + Latitude_stdz + (1|SiteName)). We then used the ‘dredge’ function in the ‘MuMIn’ R package to fit all the possible combinations of models nested in the global models. Model selection was performed using an Information-Theory approach (Burnham and Anderson, 2002), based on Akaike Information Criterion values corrected for small sample size (AICc). For each replication we fitted and ranked the global model and the submodels but only the top ranking model were chosen. We determined conditional and marginal R2 following the method of Nakagawa and Holger (2013) to estimate the explained variance of the fixed and random effects in the LMM. Conditional R2 represents the variance explained by both fixed and random effects, and marginal R2 refers to the variance explained by fixed effects only. The difference between these two components gives the R2 of the random effect. All candidate models with delta < 2 are presented (Burnham and Anderson, 2002). All analyses were performed in R 3.4.5 (R Core Team, 2018).