## 2 Materials and Methods
### Statistical Analyses  
A modified comparable soil quality index was developed to determine the effect of the road on soil quality associated with soil nutrients and soil heavy metals referring to the soil deterioration index [42; 43]. It was computed on the assumption that the status of the individual soil properties was once the same as that of adjacent soils. The equation of \(SQI\) is expressed below:  
\[SQI=\sum_{i=1}^{n}Q(x_{i})\times 100\%/n \tag{1}\]  
where \(x_{i}\) is the value of soil properties selected for the soil quality; \(Q(x_{i})\) is the membership values of different soil properties; \(Q(x_{i})\) were calculated by the ascending and descending functions. Descending function was used for heavy metals as its higher contents often indicated soil quality degradation. \[Q(x_{i}) =(x_{ij}-x_{\text{imin}})/(x_{\text{imax}}-x_{\text{imin}}) \tag{2}\] \[Q(x_{i}) =(x_{\text{imax}}-x_{ij})/(x_{\text{imax}}-x_{\text{imin}}) \tag{3}\]where \(x_{ij}\) is the values of soil properties; \(x_{imax}\) and \(x_{imin}\) are the maximum and minimum value of the \(i\) soil property in the same transect line. So, \(SQI\) indexes could be compared across all the transect lines at different locations. In addition, multivariable statistic approaches were performed for influencing factor analysis. The PCA method was used to identify soil heavy metal sources. The paired-samples \(T\) test was used to compare the metal concentrations under different land uses and near different road levels. The Kaiser-Meyer-Olkin (KMO) and Bartlett's sphericity tests were performed to examine the suitability of the data. One-way analysis of variance (ANOVA) procedures were used to compare the effects of different influencing factors. The ANOVA was conducted using the SAS program. Finally, we used CCA for direct gradient analysis to identify the relative importance of different influencing factors of heavy metals for the natural soils. It has been suggested that CCA is the most feasible ordination method for performing direct gradient analysis. The software program Canoco 4.0 was used for data processing [44]. The data sets in our study contain \(i\) plots (\(i=1\) to \(n\)), \(k\) influencing variables (\(k=1\) to \(m\)) and \(q\) soil heavy metal attributes. We selected the pH, TP and TK as the soil property indicators significantly related to heavy metals, plant layer coverage as the land cover indicator, relative elevation as the topography indicator, and distance off road as the horizontal influencing indicator. These influencing factors were set as environmental independent variables and the soil heavy metals as dependent variables. Let Y be the data matrix of the variables, an \(n\times k\) matrix of variables data (119 \(\times\) 6 matrix), E be the matrix of soil heavy metal attributes, and an \(n\times p\) matrix of dependent variables data (119 \(\times\) 7 matrix) be the plots of heavy metal attributes. Thus, \(y_{ik}\), an element of Y, is the value of variable \(k\) at plot \(i\), and \(e_{ji}\), an element of E, is the value of soil heavy metal attribute \(j\) at plot \(i\). The following steps are performed within each CCA iteration [31]:  
(1) Choose arbitrary, but unequal, initial plot scores (\(x_{i}\)). (2) Calculate soil series scores (\(u_{k}\)) by weighted averaging of the plot scores with \(\lambda=1\) (where \(\lambda\) is eigenvalue; see step 7).