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).